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Abstract 

The  finite-element  shock  hydrodynamics  code  ALEGRA  has  recently  been  upgraded  to 
include  an  X-FEM  implementation  in  2D  for  simulating  impact,  sliding,  and  release 
between  materials  in  the  Eulerian  frame.  For  validation  testing  purposes,  the  problem  of 
long-rod  penetration  in  semi-infinite  targets  is  considered  in  this  report,  at  velocities  of 
500  to  3000  m/s.  We  describe  testing  simulations  done  using  AEEGRA  with  and  without 
the  X-FEM  capability,  in  order  to  verify  its  adequacy  by  showing  X-FEM  recovers  the 
good  results  found  with  the  standard  AEEGRA  formulation.  The  X-FEM  results  for 
depth  of  penetration  differ  from  previously  measured  experimental  data  by  less  than  2%, 
and  from  the  standard  formulation  results  by  less  than  1%.  They  converge  monotonically 
under  mesh  refinement  at  first  order.  Sensitivities  to  domain  size  and  rear  boundary 
condition  are  investigated  and  shown  to  be  small.  Aside  from  some  simulation  stability 
issues,  X-FEM  is  found  to  produce  good  results  for  this  classical  impact  and  penetration 
problem. 
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NOMENCLATURE 


ALE 

Arbitrary  Eagrangian-Eulerian 

BC 

Boundary  Condition 

DOP,  P 

Depth  Of  Penetration 

EOS 

Equation  Of  State 

KEOS 

“Kerley”  Equation  Of  State 

E_ROD,  E 

Rod  (projectile)  Eength 

MB 

MegaByte 

PIR 

Pattern  Interface  Reconstruction 

R_PET 

Radius  of  target  PEaTe 

R_ROD 

Radius  of  ROD  (projectile) 

RHA 

Rolled  Homogenous  Armor 

S_CEE 

Cell  Size 

SR 

Density  scaling  factor  (“Scale  Rho”) 

SWRI 

Southwest  Research  Institute 

T_PET 

Thickness  of  target  PEaTe 

V_ROD, Vp 

Rod  (projectile)  Velocity 

X-EEM 

extended  Einite  Element  Method 
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1.  INTRODUCTION 


1.1.  Background 

The  finite-element  shock  hydrodynamics  code  ALEGRAf  [18]  has  recently  been  upgraded  to 
include  an  “extended  finite  element  method”  (X-FEM)  implementation  in  2D  for  simulating 
impact,  sliding,  and  release  between  materials  in  the  Eulerian  frame  [6,19,22,23].  This  is  a 
significant  extension  of  the  capabilities  of  AEEGRA  and  as  such  requires  testing  in 
environments  that  are  representative  of  the  software  application  space.  One  application  for 
which  AEEGRA  is  used  is  high-velocity  impact  and  penetration  events.  These  are  inherently 
multi-material  systems  which  include  strong  shock  waves  and  high  deformations.  Often,  they 
require  the  use  of  shock-capturing  methods,  material  equation  of  state  (EOS)  and  constitutive 
models,  and  Eulerian  meshes  with  remap  algorithms  that  can  maintain  the  validity  of  the 
calculation  despite  large  strains  and  strain  rates. 

The  problem  of  long-rod  penetration  in  semi-infinite  targets  is  considered  here  as  representative 
of  one  class  of  impact  and  penetration  events.  This  problem  has  been  studied  in  detail  in  the 
experiments  by  Hohler  and  Stilp  (1977,1987)  [8,9].  Their  experiments  involving  the  impact  of 
long  tungsten-heavy-alloy  cylinders  on  plates  of  hardened  steel,  reproduced  numerically  by 
Anderson,  Jr.  et  al.  (1991)  [3],  form  a  basis  for  usability  testing  and  validation  done  in  support  of 
AEEGRA  code  development  activities. 

In  the  present  report,  we  describe  testing  simulations  done  using  AEEGRA  with  and  without  the 
X-EEM  capability,  in  order  to  evaluate  the  adequacy  of  X-EEM  for  modeling  impact  and 
penetration.  After  a  brief  description  of  the  numerical  methods  used  in  this  study,  the  remainder 
of  the  report  is  organized  as  follows.  Experimental  setup  and  results  are  discussed  briefly  in 
Section  2.  The  setup  for  AEEGRA  simulations  is  described  in  detail  in  Section  3.  Results  from 
these  simulations  are  presented  in  Section  4,  with  discussion  of  post-processing  techniques, 
effects  of  the  mesh  setup  and  resolution,  and  comparison  of  results  with  and  without  X-EEM  to 
the  experimental  data. 

1.2.  Numerical  methods 

AEEGRA  is  used  for  all  numerical  simulations  in  the  present  study.  AEEGRA  is  an  arbitrary 
Eagrangian-Eulerian  shock  hydrodynamics  and  multi-physics  code  developed  and  maintained  at 
Sandia  National  Laboratories  since  the  1990’s.  AEEGRA  uses  the  finite-element  method  with 
explicit  time  integration  to  solve  the  equations  of  solid  dynamics  for  multi-material  media 
subjected  to  shocks  and  strong  deformations,  using  material-specific  equation  of  state  (EOS)  and 
constitutive  models  to  close  the  system.  AEEGRA  also  incorporates  second-order-accurate 
interface  tracking,  remap  algorithms  that  are  third-order-accurate  in  one  dimension,  modern 
artificial  viscosity  methods  for  handling  shock  waves,  and  stabilized  “isentropic”  multi-material 
treatment  with  radiative  temperature  relaxation  [17,18].  AEEGRA  validation  efforts  include  the 
work  by  Carroll  et  al.  (1997)  [5]  for  very  similar  problems  to  those  considered  here:  high- 
velocity  impact  and  penetration  of  long  rods  in  semi-infinite  targets.  Other  efforts  include  the 


t  Arbitrary  Lagrangian-Eulerian  multi-physics  finite-element  code  for  large  distortion  and  shock  propagation 
developed  by  Sandia  National  Laboratories. 
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work  by  Lemke  et  al.  (2005)  [14]  for  magnetically  accelerated  flyer  plates;  and  Doney  et  al. 
(2010)  [7]  for  exploding  wires,  among  others. 

ALEGRA  simulates  the  motion  of  material  across  the  domain  using  the  arbitrary  Lagrangian- 
Eulerian  (AEE)  approach,  which  moves  material  with  the  computational  mesh  in  an  initial 
Lagrangian  step,  then  reconfigures  the  mesh  in  a  subsequent  “remesh”  and  then  a  “remap”  step. 
Most  commonly,  the  mesh  is  reconfigured  back  to  its  original,  stationary  state,  which  is  the 
Eulerian  frame.  This  remesh  and  remap  operation  allows  simulations  to  proceed  despite  large 
shear  and  deformations,  but  requires  the  use  of  advection  algorithms  to  move  the  solution 
information  between  the  meshes. 

Since  Eulerian  simulations  necessarily  generate  elements  that  contain  multiple  materials, 
additional  infrastructure  is  required  to  reconstruct  the  location  and  orientation  of  material 
interfaces,  to  partition  the  deformations  to  individual  materials  in  an  element,  and  to  construct  a 
total  element  stress  based  on  stresses  computed  for  individual  materials  in  an  element.  All  of 
these  operations  introduce  significant  approximations  and  errors  into  the  simulations,  resulting  in 
spurious  behavior  in  some  situations.  Eor  example,  since  the  method  only  allows  a  single  stress 
and  deformation  in  an  element,  two  materials  within  an  element  cannot  separate  and  allow  a  gap 
to  open  when  subjected  to  tensile  stresses.  Thus  the  materials  are  effectively  bonded,  and  only 
the  use  of  ad  hoc  modifications  to  the  material  state  can  allow  them  to  separate. 

The  X-EEM  has  been  implemented  in  AEEGRA  in  2D  in  order  to  eliminate  these  issues,  so  that 
inter-material  dynamics  can  be  handled  as  surface  phenomena  even  in  the  Eulerian  frame, 
without  volumetric  modifications  to  the  thermodynamic  state  of  the  materials.  In  X-EEM,  the 
finite-element  basis  functions  in  each  multi-material  Eulerian  element  are  “enriched”  using 
Heaviside  functions,  so  that  individual  materials  have  basis  functions  supported  only  within  the 
region  where  the  material  exists.  Each  material  within  an  element  can  then  have  an  independent 
displacement  field,  with  interfaces  between  them  represented  explicitly  in  the  solution.  With  this 
numerical  infrastructure,  materials  can  slide  with  respect  to  or  separate  from  one  another  inside 
an  Eulerian  element.  The  concept  of  Heaviside  enrichment  is  sketched  notionally  in  Eigure  1 . 


Figure  1 :  Heaviside  enrichment  in  X-FEM.  Nodes  are  indexed  by  J,  and  materials  by  m, 
within  a  single  element  in  the  computational  mesh.  Interface  intersections  are  green. 

With  the  Heaviside  enrichment  in  place,  the  X-FEM  implementation  in  AEEGRA  employs  three 
methods  to  capture  material  motion  across  an  Eulerian  mesh.  First,  “pattern  interface 
reconstruction”  (PIR)  is  used  to  capture  material  interfaces  with  maximum  fidelity.  This  method 
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reproduces  linear  interfaces  in  all  orientations  and  has  second-order  accuracy  for  planar  and 
curved  interfaces  (when  interface  “smoothing”  is  applied).  Second,  forces  on  enriched-element 
nodes  are  controlled  by  a  forward-increment  Lagrange  multiplier  contact  algorithm  which 
enforces  interface  physics.  This  allows  inter-material  physics  to  be  treated  as  a  true  surface 
phenomenon.  Third,  remap  is  performed  using  a  monotonicity-preserving  “intersection  remap” 
scheme  with  third-order  limiting,  which  conservatively  remaps  materials  and  material  interface 
intersections  from  the  deformed  Lagrangian  mesh  to  the  Eulerian  mesh. 

Although  third-order  accuracy  in  intersection  remap  is  available,  only  the  first-order  scheme  has 
been  fully  tested  for  use  in  applications,  so  this  scheme  is  used  all  simulations  shown  here.  In 
the  standard,  “unenriched”  ALEGRA  formulation,  third-order  limiting  is  used  by  default  in 
single-material  regions,  but  remap  uses  a  “swept- volume”  approach,  rather  than  the  “intersection” 
scheme  used  with  X-EEM.  In  the  unenriched  formulation  only,  the  order  of  accuracy  used  in 
remap  can  be  controlled  directly  by  the  user,  and  is  reduced  adaptively  by  the  code  when  multi¬ 
material  elements  are  present  in  the  local  stencil.  In  the  donor  remap  configuration,  the  code 
uses  a  first-order-accurate  swept-volume  scheme  for  all  stencils.  In  the  aggressive  remap 
configuration,  the  code  uses  up  to  third-order  accuracy,  depending  on  the  stencil.  Both  of  these 
schemes  will  be  used  only  in  unenriched  (non-X-EEM)  simulations  shown  here. 

With  X-EEM,  AEEGRA  is  able  to  simulate  multi-material  shock  hydrodynamics  in  the  Eulerian 
frame,  while  capturing  the  details  of  material  interface  interactions  locally.  It  is  anticipated  that 
these  methods  will  dramatically  enhance  the  fidelity  with  which  certain  shock-hydrodynamic 
phenomena  can  be  simulated.  For  the  impact/penetration  scenario  considered  here,  ALEGRA’s 
capability  using  the  standard  formulation  is  well  established  [5].  Thus,  our  intent  is  not  to 
characterize  any  enhancement  of  fidelity,  but  to  investigate  whether  AEEGRA  still  captures  the 
details  of  classical  problems  correctly  with  X-EEM.  That  is,  we  aim  to  show  that  with  X-EEM, 
results  computed  using  the  standard  “unenriched”  Eulerian  methods  can  still  be  recovered.  To  do 
this,  we  examine  a  classical  long-rod  penetration  problem. 


2.  EXPERIMENTAL  DATA 

The  penetration  mechanism  of  long  steel-  and  high-density  rods  into  semi-infinite  steel  targets  of 
different  ultimate  yield  strength  has  been  investigated  experimentally  in  previous  work  by 
Hohler  and  Stilp  (1977,  1987)  [8,9]  at  impact  velocities  of  approximately  500  to  3000  m/s.  In  the 
component  of  Hohler  and  Ship’s  work  considered  here,  all  rods  were  tungsten  alloy  and  circular 
in  cross-section  with  a  length-to-diameter  ratio  of  LID  =  10.  The  calibers  varied  between  D  =  2.5 
to  6  mm,  the  masses  between  m  =  1  to  31  g.  The  tests  were  carried  out  by  means  of  a  powder  gun 
and  a  two-stage  light-gas  gun.  Eigure  2  shows  the  measured  penetration  depth  P  normalized  by 
the  projectile  length  L  in  the  “semi-infinite”  target  as  a  function  of  the  projectile  velocity  Vp 
obtained  from  the  experiment  [8,2]. 

These  data  appear  in  Southwest  Research  Institute  (SWRI)  Penetration  Mechanics  Database  [2] 
in  Eigure  2.1.25,  and  in  the  associated  data  tables  in  Appendix  A  (pages  A- 119  to  A- 124,  from 
the  "Hohler  and  Stilp  (1991a)"  section).  The  data  used  here  were  extracted  for  the  LID  =10  case, 
where  right  circular  cylinders  were  shot  at  normal  incidence.  The  data  do  not  include  any 
information  on  measurement  uncertainty  or  other  potential  sources  of  error,  and  no  attempt  was 
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made  here  to  quantify  the  errors.  An  inherent  variability  of  a  few  percent  is  apparent.  These 
experiments  have  been  simulated  previously  using  many  other  codes,  with  comparison  to 
experimental  data,  as  shown  by  Anderson,  Jr.  et  al.  (1991)  [3]  and  Templeton  et  al.  (2002)  [21], 
among  others.  The  experiments  have  not  been  simulated  previously  using  any  X-FEM 
implementation. 

The  projectile  material  is  the  “D17.6”  tungsten  alloy:  ambient  density  17,600  kg/m  and  Brinell 
hardness  number  406  [2].  The  target  material  is  the  “HzB,A”  armor  steel:  ambient  density  7,850 
kg/m^  and  Brinell  hardness  number  255  [2].  These  results  and  others  from  experiments  very 
similar  to  these  appear  -  in  graphical  form  only  -  in  Hohler  and  Stilp  (1977,  1987)  [8,9].  The 
D17.6  tungsten  alloy  was  used  also  in  the  time -resolved  penetration  experiments  of  Anderson  et 
al.  (1995)  [1],  who  describe  it  as  a  sintered  alloy  with  the  composition  92.5%  W,  4.85%  Ni,  2.4% 
Fe,  0.25%Co.  The  HzB,A  material  is  a  high-hard  armor  steel,  but  the  details  of  its  composition 
are  unknown.  In  the  original  experiments,  it  is  not  clear  how  the  penetration  depth  was  measured, 
but  unpublished  documents  from  Ernst  Mach  Institute  indicate  that  this  type  of  experiment 
typically  used  post-test  sectioning  of  the  target  plates  and  direct  measurement  of  the  penetration 
channel  depth  relative  to  the  original  surface  plane.  (The  later  experiments  of  Anderson  et  al. 
(1995)  [1]  used  time-resolved  x-ray  imaging.)  Hohler  and  Stilp  (1977,  1987)  note  that  the  “S- 
shape”  of  the  penetration  trend  versus  impact  velocity  is  a  consequence  of  the  hydrodynamic 
character  of  the  penetration  event  at  higher  impact  velocities  [8,9],  and  indicates  the  “dynamic 
strength”  of  the  target  material  exceeds  that  of  the  projectile  material. 


Projectile  striking  velocity,  Vp  (m/s) 

Figure  2:  Experimental  data  from  a  penetration  mechanics  database  [2]:  the  penetration 
depth  P  normalized  by  the  projectile  length  L  as  a  function  of  the  projectile  velocity  Vp, 
with  empirical  data  fits  based  on  the  Lanz-Odermatt  formula  [13]  and  a  5th-order 
polynomial. 
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A  semi-analytical  empirical  fit  to  Hohler  and  Stilp’s  data  is  included  in  Figure  2,  based  on  the 
commonly  used  inverse-exponential  form  originally  devised  by  Lanz  and  Odermatt  (1992)  [13] 
for  long-rod  penetration  in  semi-infinite  targets,  as  modified  by  Rapacki  et  al.  (1995)  [15]: 


(1) 


In  this  formula,  the  target  material  penetration  resistance  is  represented  by  the  quantity 


5  =  k{BHN)^ 


(2) 


which  has  units  of  stress.  The  projectile  and  target  ambient  densities  are  denoted  by  pp  and  p-p, 
respectively,  and  the  Brinnell  hardness  of  the  target  material  by  BHN.  In  these  formulae,  the 
quantities  a,  k,  and  M  are  obtained  by  a  least-squares  fit  of  the  formulae  to  the  experimental  data. 
The  asymptotic  normalized  penetration  at  very  high  velocities  is  represented  by  the 
“hydrodynamic  multiplier,”  a.  The  exponent  in  the  formula  represents  the  ratio  of  material 
penetration  resistance  to  projectile  dynamic  pressure.  (See  also  Reference  20.) 

For  the  “Lanz-Odermatt”  fit  shown  in  Figure  2  and  elsewhere  in  this  report,  the  target  hardness  is 
assumed  to  be  BHN  =  295,  which  is  the  mean  value  in  the  range  quoted  by  Hohler  and  Stilp 
(1977)  [8],  rather  than  the  value  quoted  in  the  Penetration  Mechanics  Database  [2].  Least- 
squares  fitting  to  the  data  yields  the  following  fit  parameters:  a  =  1.22817,  k  =  1.63084  GPa,  and 
M  =  0.373287.  These  values  compare  favorably  with  the  values  obtained  by  Rapacki  et  al.  (1995) 
[15]  for  similar  materials  in  a  separate  set  of  experiments.  However,  the  fit  fails  to  capture  data 
points  accurately  at  the  extrema  of  the  velocity  range,  and  slightly  diminishes  the  slope  of  the 
curve.  Therefore,  a  purely  empirical  fit  is  also  included  using  a  simple  5*-order  polynomial. 


3.  MODEL  DESCRIPTION 

For  validation  testing  of  X-FEM  in  ALEGRA,  simulations  for  the  test  cases  in  Table  1  were 
conducted.  The  depths  of  penetration  (DOP  =  P),  with  six  rod  velocities  (V_ROD  =  Vp),  for  each 
case  were  calculated  as  a  result  from  an  ALEGRA  run.  Each  DOP  was  normalized  by  the 
original  length  of  the  projectile  (L_ROD  =  L).  The  normalized  penetration  depths  (PIL)  with 
V_RODs  for  the  baseline  model  were  compared  to  the  experimental  data. 

Eigure  3  shows  the  schematic  diagram  of  the  model  and  boundary  conditions  (BC).  The  target 
plate  thickness  (T_PLT)  and  radius  (R_PLT)  for  the  baseline  are  120  mm  and  40  mm, 
respectively.  The  projectile  rod  length  (L_ROD)  and  radius  (R_ROD)  for  the  baseline  are  50  mm 
and  2.5  mm,  respectively.  To  investigate  the  model  parameter  effect,  the  boundary  condition  on 
the  plate  back  surface,  T_PLT,  and  R_PLT  were  changed  as  listed  in  Table  1,  using  either  a  free- 
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motion  (unconstrained)  condition,  or  a  “no-displacement”  condition,  which  prevents  motion 
normal  to  the  boundary.  The  denotation  of  each  parameter  in  Table  1  is  shown  in  Figure  3. 


Table  1 :  The  parameter  values  for  each  test  case. 


Test  case 

Rod  radius, 

R  ROD  (mm) 

Rod  length, 

L  ROD  (mm) 

Plate  radius, 
R PLT  (mm) 

Plate  thickness, 
L PLT  (mm) 

BC  on  Plate 
Back  Surface 

4x12fixed 

2.5 

50. 

40. 

120. 

No  displacem’t 

4x12free 

2.5 

50. 

40. 

120. 

Free 

4x18fixed 

2.5 

50. 

40. 

180. 

No  displacem’t 

5x12fixed 

2.5 

50. 

50. 

120. 

No  displacem’t 

Note:  4x12fixed  is  the  baseline. 


Center 

Line 


Figure  3:  Schematic  diagram  and  boundary  conditions  of  the  model  with  denotation  of 
each  parameter. 
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Making  use  of  axisymmetry,  only  a  half-section  of  the  target  plate  and  the  rod  projectile  are 
modeled  as  shown  Figure  3.  The  projectile  moves  from  bottom  to  top,  in  the  -i-z-direction.  The 
left  boundary  is  the  axis  of  axisymmetry.  A  zero-displacement  boundary  condition  (radial 
direction)  is  applied  on  the  right  boundary  of  the  model  to  represent  the  “semi-infinite”  steel 
target.  A  zero-displacement  boundary  condition  (axial  direction)  is  applied  on  the  upper 
boundary  to  better  represent  the  projectile  impact  conditions  in  the  experiments  for  the  baseline 
model. 

The  BC’s  and  dimensions  of  the  target  were  chosen  to  represent  impact  and  penetration  in  a 
semi-infinite  target.  The  actual  dimensions  of  the  experimental  target  are  unknown,  but  all  of  the 
data  are  categorized  in  the  Penetration  Mechanics  Database  under  “Semi-Infinite  Data”  [2], 
meaning  that  the  target  plate  was  large  enough  to  stop  all  projectile  impacts  for  the  given 
velocities,  without  spall  or  breakout  on  the  rear  surface.  Since  we  cannot  model  the  entire 
dimensions  of  the  real  plate,  we  must  model  the  semi-infinite  effect  as  best  we  can  using  a  plate 
of  finite  dimensions  and  fixed  or  free  boundary  conditions.  The  purpose  of  the  runs  listed  in 
Table  1  was  to  identify  which  set  of  target  dimensions  and  BC’s  most  accurately  represents  the 
semi-infinite  system. 

The  equation  of  state  (EOS)  material  properties  for  the  projectile  rod  and  target  plate  were 
obtained  from  Sesame  tabular  libraries,  using  the  “KEOS”  (Kerley)  interface  in  AEEGRA,  as 
listed  in  Table  2.  The  KEOS  Sesame  model  is  an  implementation  of  the  Sesame  EOS  in 
AEEGRA.  The  KEOS  Sesame  model  uses  the  same  tables  that  are  used  by  the  CTHJ  SES  model. 
The  von  Mises  flow  stress  model  is  used  for  the  yield  calculation.  The  model  for  the  von  Mises 
flow  stress,  a,  is  the  form  of  Johnson  and  Cook  (1983)  [10]: 


a  =  [A  +  Be^][l  +  C\n£*][l-T*^]  (3) 

where  e  is  the  equivalent  plastic  strain,  e*  —  e/ e'g  is  the  dimensionless  plastic  strain  rate  for 
do  =  1.0  s“^,  and  T*  is  the  homologous  temperature,  which  quantifies  the  departure  from  the 
melt  temperature  Tm.  The  quantities  A,  B,  n,  C,  and  m  are  material  constants  obtained 
empirically  for  common  materials.  The  expression  in  the  first  set  of  brackets  gives  the  stress  as  a 
function  of  strain  for  e*  =  1.0  and  T*  —  0.  The  expressions  in  the  second  and  third  sets  of 
brackets  represent  the  effects  of  strain  rate  and  temperature,  respectively.  The  Johnson-Cook 
model  parameters  used  in  AEEGRA  in  this  study  are  shown  in  Table  2. 

The  yield  model  and  parameters  used  here  are  those  used  in  CTH  and  EPIC  by  Templeton  et  al. 
(2002)  [21].  As  these  parameters  were  used  successfully  in  this  previous  work  (<10%  deviation 
from  Hohler  and  Stilp’s  data),  they  were  used  again  here.  However,  no  conclusive  evidence 
exists  which  establishes  any  certain  combination  of  material  response  models  for  this  problem  as 
the  most  accurate.  The  tabular  Sesame  EOS  models  are  used  here,  rather  than  the  analytic  Mie- 
Griineisen  EOS  models  used  previously.  Other  successful  parameter  settings  for  this  problem  or 
similar  problems  and  materials  can  be  found  in  Kmetyk  and  Yarrington  (1988)  [12];  Anderson, 
Jr.  and  Walker  (1991)  [3];  Anderson,  Jr.  et  al.  (1995)  [1];  and  Carroll  et  al.  (1997)  [5]. 


J  CTH  is  a  multi-material,  large  deformation,  strong  shock  wave,  solid  mechanics  code  developed  at  Sandia 
National  Laboratories. 
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Table  2:  Material  properties  used  in 


he  analysis. 


Input 

Parameter 

Unit 

Projectile  Rod 
(Tungsten) 

Target  Plate 
(HzB,  Steel) 

KEOS  Model 

FEOS 

- 

SESAME 

SESAME 

NEOS 

- 

3550 

2150 

Density  (p) 

kg/m^ 

17760 

7850 

SR 

- 

1.094 

1.003 

Yield  Model 
(Johnson- 
Cook) 

A 

GPa 

1.365 

0.810 

B 

GPa 

0.1765 

0.5095 

C 

- 

0.016 

0.014 

m 

1.00 

1.030 

n 

0.12 

0.260 

Tm 

K 

3695 

1818 

Poisson  (v) 

- 

0.281 

0.299 

Note:  -  KEOS  stands  for  “Kerley”  EOS.  KEOS  is  an  interface  to  Sesame  data  that  was  built  by  Gerald  Kerley  and  is 
very  commonly  used  in  ALEGRA  simulations. 

-  SR  is  a  density  scaling  factor,  which  represents  the  ratio  of  molecular  weight  for  the  table  relative  to  the 
actual  material.  For  alloys  it  is  commonly  used  to  relate  the  pure  material  density  to  the  alloy  density. 


4.  ANALYSIS  RESULTS 
4.1 .  Depth  of  penetration 

To  calculate  the  depth  of  penetration  (DOP)  from  ALEGRA  simulations,  Lagrangian  “Tracers” 
H  and  P  were  placed  at  the  center  of  the  rod  nose  or  “head”  and  at  the  edge  of  the  target  or  “plate” 
front  surface,  respectively,  as  shown  in  Figure  3.  “Tracers”  are  diagnostic  markers  used  to  track 
material  or  mesh  motion  locally  in  an  ALEGRA  calculation.  Lagrangian  tracers  are  advected 
with  the  medium  by  interpolation  of  the  velocity  field.  Figure  4  shows  the  z-positions  of  Tracers 
H  (“rod  head”)  and  P  (“plate  surface”)  as  a  function  of  time  from  the  baseline  model  with 
V_ROD=3000  m/s  as  an  example.  Figure  5  presents  views  of  progression  of  the  formation  of  the 
penetration  channel  at  0  and  49  ps.  We  observe  that  the  rod  nose/head  position  goes  through  a 
maximum,  and  subsequently  recoils.  Also,  the  plate  surface  near  the  penetration  channel  moves 
downward  (-z-direction)  as  the  rod  head  moves  into  the  target,  due  to  crater  formation.  The  plate 
surface  position  at  the  outer  radius  continues  to  oscillate  through  the  duration  of  the  simulation. 

To  mimic  a  post-test  penetration  measurement,  DOP  data  extracted  from  the  ALEGRA 
simulations  must  use  the  maximum  z-position  reached  by  the  rod  head,  and  some  representation 
of  the  final  z-position  of  the  plate  surface,  after  oscillations  have  damped  out.  To  approximate 
this,  we  take  the  perpendicular  distance  between  the  positions  of  rod  head  (Tracer  H)  and  the 
plate  front  surface  (Tracer  P)  at  the  time  of  maximum  Tracer-H  z-position  to  be  the  DOP. 
(Continuing  the  simulation  until  plate  oscillations  damp  out  would  be  prohibitively  expensive.) 
In  the  V_ROD=3000  m/s  case  shown  in  Figure  4,  the  positions  of  rod  head  and  plate  surface  are 
calculated  to  be  74.5  mm  and  -1.8  mm,  respectively,  at  49  ps  (time  at  the  peak  of  the  Tracer  H 
curve).  Thus,  the  DOP  is  calculated  to  be  76.3  mm. 
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Tracer  z-position  (m) 


Figure  4:  Calculation  of  OOP  for  the  baseline  case  with  V_ROD=3000  m/s. 
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Figure  5:  Simulation  results  at  0  and  49  ps  for  the  baseline  case  with  V_ROD=3000  m/s 
(Unit  on  axes  is  m). 
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4.2.  Mesh  resolution 


Use  of  the  X-FEM  algorithm  in  ALEGRA  instead  of  the  standard  Eulerian  formulation  incurs 
dramatically  increased  computational  cost.  The  mesh  with  a  smaller  cell  size  provides  higher 
accuracy,  but  spends  more  computer  running  time.  Therefore,  an  appropriate  mesh  cell  size  and 
configuration  must  be  determined.  To  examine  the  cell  size  effect  on  the  accuracy  of  the  analysis 
result  and  the  computer  running  time,  meshes  with  various  cell  size  are  constructed,  shown  in 
Eigure  6  through  Eigure  11.  The  DOPs  are  calculated  from  the  AEEGRA  solutions  using  these 
meshes,  with  and  without  X-EEM.  The  geometries  of  the  baseline  (4xl2fixed)  in  Table  1  and 
material  properties  in  Table  2  are  used  as  input  parameter  values.  A  projectile  velocity  of  1500 
m/s  is  used  for  all  cases  for  comparison  purposes.  The  AEEGRA  runs  are  performed  using  64 
nodes  (512  processors)  on  RedSky^ 

The  meshes  in  Eigure  6  through  Eigure  8  were  generated  using  AEEGRA’ s  inline  meshing 
capability,  as  uniform  meshes.  The  cell  sizes  in  the  x-  and  y-  directions  are  the  same  as  the  input 
cell  size  value.  The  projectile  rod  head  is  expected  to  move  between  0  mm  to  78  mm  on  the 
center  line  of  the  model  since  the  penetration  depth  was  measured  as  74.1  mm  for  the  projectile 
velocity  of  3000  m/s  (maximum  velocity  in  this  study)  [2].  The  crater  radius  is  expected  to  be 
less  than  four  times  the  projectile  rod  radius.  The  region  0  <  z  <  78  mm  and  0  <  r  <  10  mm  is 
defined  as  the  “penetration  area”.  The  same  size  of  cells  in  the  outside  and  inside  of  the 
penetration  area  is  inefficient  for  the  computer  run.  The  graded  cell  size  meshes  as  shown  Eigure 
9  through  Eigure  11  were  generated  using  the  programmable  cell  size  options  in  AEEGRA’ s 
inline  mesh  capability.  Cells  smaller  than  the  input  cell  size  value  (S_CEE)  were  generated  in 
the  penetration  area.  In  the  5-G  subcase,  the  cell  size  input  value  of  0.15  mm  generated  the  cell 
sizes  of  0.05  mm  and  0.07  mm  in  the  r-  and  z-  directions.  The  cell  size  gradually  increases 
outside  the  penetration  area.  The  number  of  elements  in  5-G  mesh  is  smaller  than  that  of  3-U 
mesh  while  the  cell  size  for  the  5-G  subcase  in  the  penetration  area  is  smaller  than  that  of 
subcase  3-U:  i.e.,  subcase  5-G  will  produce  more  accurate  results  with  fewer  elements  than 
subcase  3-U,  due  to  the  graded  mesh. 


§  Red  Sky  is  a  platform  which  is  assembled  in  the  space  where  the  legendary  system  ASCI  Red  once  stood  at  Sandia 
National  Laboratories. 
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Figure  6:  Computational  mesh  for  uniform  cell  size=0.5  mm  (Unit  on  axes:  m). 
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Figure  7:  Computational  mesh  for  uniform  cell  size=0.25  mm  (Unit  on  axes:  m). 
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Figure  8:  Computational  mesh  for  uniform  cell  size=0.125  mm  (Unit  on  axes:  m). 
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Figure  9:  Computational  mesh  for  graded  cell  size=0.2  mm  (Unit  on  axes:  m). 
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Figure  10:  Computational  mesh  for  graded  cell  size=0.15  mm  (Unit  on  axes:  m). 
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Figure  11 :  Computational  mesh  for  graded  cell  size=0.1  mm  (Unit  on  axes:  m). 


21 


A  series  of  simulations  was  performed  at  varying  mesh  resolution  for  1500-m/s  impact,  using 
both  ALEGRA’s  X-FEM  and  standard  “unenriched,”  or  “U-EEM”  methods.  Uniform  and 
graded  meshes  were  used  to  achieve  the  desired  mesh  resolution  in  the  penetration  region.  Table 
3  lists  the  mesh  type,  cell  size,  number  of  elements,  size  of  EXODUS  II**  file  as  an  output, 
actual  elapsed  time  for  each  of  the  X-EEM  computer  runs,  and  calculated  DOP  for  each  mesh. 


Table  3:  X-FEM  simulation  results  at  1500  m/s  for  varying  mesh  resolutions. 


Subcase 

Type 

Input  cell 
size, 

S  CEL 
(mm) 

Cell  size  in 
penetration  area 
(mm) 

No.  of 
elements 

EXODUS 
File  Size 
(MB) 

Elapsed 

time 

(hh:mm: 

ss) 

Predicted 

DOP 

(mm) 

Figure 

Ar 

Az 

h 

1-U 

Uniform 

0.5 

0.5 

0.5 

0.5 

27280 

375 

00:27:21 

38.134 

Figure  6 

2-U 

Uniform 

0.25 

0.25 

0.25 

0.25 

109280 

1500 

00:44:54 

41.895 

Figure  7 

3-U 

Uniform 

0.125 

0.125 

0.125 

0.125 

437440 

6028 

02:23:38 

43.696 

Figure  8 

4-G 

Graded 

0.2 

0.067 

0.129 

0.098 

170800 

2355 

02:15:36 

44.004 

Figure  9 

5-G 

Graded 

0.15 

0.052 

0.073 

0.063 

301644 

3316 

04:09:38 

44.546 

Figure  10 

6-G 

Graded 

0.1 

0.035 

0.049 

0.042 

680000 

5852 

08:20:08 

44.825 
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Figure  12:  Predicted  DOP,  rod  head  and  plate  surface  histories  with  VP=1500  m/s  from  X- 
FEM  simulation  using  the  mesh  with  various  cell  sizes. 


**  EXODUSII  is  a  file  format  developed  to  store  and  retrieve  data  for  finite  element  analyses.  It  is  used  for  pre¬ 
processing  (problem  definition),  post-processing  (results  visualization),  as  well  as  code-to-code  data  transfer.  An 
EXODUS  II  data  file  is  a  random  access,  machine  independent,  binary  file  that  is  written  and  read  via  C,  C-H-,  or 
EORTRAN  library  routines  which  comprise  the  Application  Programming  Interface  (API). 
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As  described  above,  transient  penetration  data  for  each  of  these  simulations  were  obtained  using 
curves  of  rod  nose/head  (Tracer  H)  and  plate  surface  (Tracer  P)  for  the  subcases.  These  histories 
are  shown  graphically  in  Figure  12.  The  data  indicate  that  as  the  mesh  is  refined,  the  DOP  and 
time  of  DOP  both  increase,  but  are  convergent.  The  motion  of  the  front  plate  contributes 
relatively  little  to  the  DOP. 

These  data  are  rearranged  and  overlaid  on  the  experimental  data  and  fits  in  Figure  13,  which 
shows  the  calculated  normalized  penetration  depths  {PIL)  obtained  using  ALEGRA  with  X-FEM 
for  the  projectile  velocity  of  1500  m/s  from  the  meshes  with  various  cell  sizes.  Table  3  lists  the 
mesh  type,  cell  size,  number  of  elements,  size  of  EXODUS  II  file  as  an  output,  actual  elapsed 
time  for  each  of  the  X-EEM  computer  runs,  and  calculated  DOP  for  each  mesh.  The 
characteristic  element  dimension  shown  in  Table  3  for  uniform  and  graded  meshes  is  given 
by  /i  =  (  Ar  +  Az)/2. 


0  500  1000  1500  2000  2500  3000  3500 

Projectile  striking  velocity,  Vp  (m/s) 

Figure  13:  Depths  of  penetration  with  Vp=1500  m/s  from  X-FEM  simulation  using  the 
mesh  with  various  cell  sizes. 

We  see  qualitatively  that  the  DOP  increases  toward  the  experimental  results  as  the  mesh  is 
refined.  To  quantify  this  convergence  trend,  the  following  technique  is  used,  based  on  the 
“multi-regression”  technique  of  Rider  and  Kamm  (2012)  [16]  and  Kamm  et  al.  (2012)  [11].  We 
begin  by  postulating  that  the  simulated  DOP  data  can  be  represented  in  the  form  P  =  a  +  bh^, 
which  appears  as  Equation  15  in  Reference  16.  Eeast-squares  fitting  of  this  expression  to  the 
simulation  results  is  used  to  infer  both  an  intrinsic  convergence  rate  and  an  asymptotic 
“converged”  value  of  P.  The  rate  appears  as  c  in  this  expression,  and  the  converged  value  is  a. 
The  least-squares  fit  produces  estimates  for  these  values,  given  the  simulation  DOP  data. 
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Uncertainty  bounds  on  the  convergence  rate  are  obtained  by  assuming  an  error  of  +h  in  the  DOP 
value.  We  define  upper  and  lower  bounds  on  the  simulated  P  values  as  P+  =  P  +  h  and 
P_  =  P  —  h,  and  then  repeat  the  least-squares  fit  for  P+and  P_  each  as  a  function  of  h.  This 
provides  uncertainty  bounds  on  the  inferred  convergence  trend. 

This  technique  yields  the  three  curves  shown  in  Figure  14.  Since  the  convergence  trend  is 
smooth  and  monotonic,  the  “screening”  recommended  by  Rider  and  Kamm  (2012)  [16]  is  not 
needed.  The  inferred  mean  convergence  rate  is  c  =  1.061,  with  an  uncertainty  of  -I-0.005/-0.004. 
The  X-FEM  results  thus  demonstrate  extremely  tight  adherence  to  a  smooth,  monotonic,  first- 
order  convergence  trend,  as  expected  for  this  class  of  problems.  Furthermore,  the  “converged” 
value  of  P  produced  by  this  analysis  is  a  =  45.33  mm.  Evaluating  the  polynomial  fit  shown  in 
Figure  2  at  Vp=1500  m/s  yields  a  reference  experimental  DOP  of  45.42  mm.  The  converged  X- 
FEM  simulation  result  therefore  differs  from  the  experimental  expectation  by  only  0.2%.  Very 
similar  convergence  behavior  was  obtained  for  U-EEM  results  as  well  (not  shown  here). 


Figure  14:  Mesh  sensitivity  for  X-FEM  and  standard  ALEGRA  ("U-FEM")  results. 


It  should  be  noted  that  the  simulation  for  an  input  cell  size  of  0.1  mm  (6-G)  ended  after  125  ps, 
as  can  be  seen  by  examination  of  Figure  12.  Failures  such  as  this  persist  but  only  affect  a 
relatively  small  number  of  X-FEM  simulations  (see  Sections  4.4. 1-4.4. 3).  Since  the  X-FEM 
implementation  is  very  recently  developed,  it  is  anticipated  that  ongoing  and  future  work  will 
improve  the  robustness,  stability,  and  performance  of  AEEGRA  simulations  with  X-FEM. 

Due  to  the  large  costs  and  risk  of  simulation  stopping  involved  with  the  finest  mesh  (6-G),  the 
remainder  of  this  study  was  conducted  using  the  graded  cell  size  of  0.15  mm  as  the  most 
appropriate  and  expedient  mesh  configuration  (5-G,  average  cell  dimension  in  penetration  area:  h 
=  0.063  mm).  This  is  the  only  configuration  which  allows  runtimes  under  5  hours  and  DOP  at 
1500  m/s  within  2%  of  the  converged  value. 
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4.3.  Comparison  to  the  standard  finite-element  algorithms 


For  validation  testing  purposes,  DOPs  are  calculated  using  X-FEM  simulations  for  projectile 
velocities  of  500,  1000,  1500,  2000,  2500,  and  3000  m/s,  with  the  baseline  model  (4xl2fixed) 
and  the  5-G  subcase.  Figure  15  shows  the  Tracer  H  (rod  head)  and  Tracer  P  (plate  surface) 
histories  for  six  velocities  of  projectile.  The  DOPs  are  predicted  to  be  1.98,  17.18,  44.55,  63.47, 
71.95,  and  76.31  mm  for  the  projectile  velocities  of  500,  1000,  1500,  2000,  2500,  and  3000  m/s, 
respectively.  The  predicted  DOPs  normalized  by  the  projectile  length  are  very  close  to  the 
experimental  data  as  shown  Figure  16. 

For  comparison  purposes,  the  DOPs  are  calculated  using  aggressive  and  donor  advection  in  the 
standard  ALEGRA  “unenriched”  finite  element  formulation.  Henceforward,  the  standard 
unemiched  formulation  will  be  referred  to  as  “U-FEM.”  Table  4  and  Eigure  17  summarize  the 
results  from  three  simulations  with  the  polynomial  fit  to  the  experimental  data.  The  X-FEM 
algorithm  yields  the  closest  DOP  to  the  data  when  the  projectile  velocity  is  3000  m/s,  while  the 
U-EEM  aggressive  algorithm  yields  the  closest  DOP  for  projectile  velocities  of  1000,  1500,  2000, 
and  2500  m/s;  and  the  donor  algorithm  yields  the  closest  DOP  for  projectile  velocity  of  500  m/s. 
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Figure  15:  Predicted  DOP,  rod  head  and  plate  surface  histories  from  X-FEM  with  six 
velocities  of  rod. 
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Figure  16:  Comparison  between  the  experimental  data  with  predicted  normalized 
penetration  depths  from  X-FEM  simulation. 
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Figure  17:  Comparison  between  the  predicted  normalized  penetration  depths  from  X-FEM, 
aggressive  and  donor  aggressive  simulations,  and  experimentai  data. 
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In  general,  all  three  algorithms  compare  very  favorably  to  the  experimental  data.  When 
computed  in  terms  of  mean  difference  in  PIL  for  all  six  velocities,  the  X-FEM  data  differ  from 
the  polynomial  fit  to  experimental  data  by  1.58%.  They  differ  from  the  U-FEM  aggressive 
results  by  0.51%.  The  U-FEM  aggressive  results  differ  by  1.40%  from  the  experimental  PIL 
polynomial  fit,  and  the  U-FEM  donor  results  by  1.62%.  These  differences  (shown  in  Table  4) 
are  well  within  any  plausible  measurement  uncertainty  that  may  be  assumed  to  be  inherent  to  the 
experimental  data.  It  is  interesting  to  note  that  the  penetration  depths  increase  toward  the 
experimental  data  as  the  remap  method  in  U-FEM  is  changed  from  donor  (first-order)  to 
aggressive  (third-order).  The  “intersection  remap”  algorithm  used  in  this  study  for  the  X-FEM 
simulations  in  this  study  is  only  first-order,  as  noted  in  Section  1.2. 

Table  4:  Predicted  DOP  and  P/L  with  six  projectile  velocities  with  polynomial-fit 
experimental  data.  (The  bold  fonts  indicate  the  closest  result  to  the  experimental  data  for 
each  velocity). 


Projectile 

X-FEM 

UFEM, 

aggressive 

UFEM,  donor 

Data 

(polynomial  fit) 

Velocity 

(m/s) 

DOP 

(mm) 

P/L 

DOP 

(mm) 

P/L 

DOP 

(mm) 

P/L 

DOP 

(mm) 

P/L 

500 

1.98 

0.040 

1.96 

0.039 

1.98 

0.040 

2.07 

0.041 

1000 

17.18 

0.344 

17.52 

0.350 

17.10 

0.342 

18.55 

0.371 

1500 

44.55 

0.891 

45.48 

0.910 

44.73 

0.895 

45.42 

0.908 

2000 

63.47 

1.269 

63.90 

1.278 

63.48 

1.270 

63.82 

1.276 

2500 

71.95 

1.439 

72.08 

1.442 

71.87 

1.437 

72.06 

1.441 

3000 

76.31 

1.526 

76.53 

1.531 

76.44 

1.529 

73.75 

1.475 

Mean 

difference 

1 .58% 

1 .40% 

1 .62% 

Note:  UFEM  =  ’’Unenriched”  FEM  =  Standard  ALEGRA. 


Table  5  lists  the  actual  computer  run  elapsed  time  for  each  simulation  with  the  projectile 
velocities.  The  AEEGRA  runs  were  performed  using  64  nodes  (512  processors)  on  RedSky.  The 
X-FEM  algorithm  spent  longer  time  than  the  standard  U-FEM  algorithms  by  a  factor  of  3  to  4. 
The  donor  algorithm  spent  relatively  shorter  time  than  others,  with  similar  accuracy.  (This 
outcome  is  not  expected  for  all  classes  of  problems.) 

Table  5:  Actual  elapsed  time  for  the  computer  run. 


Projectile 

Velocity 

(m/s) 

X-FEM 

(hh:mm:ss) 

UFEM, 

aggressive 

(hh:mm:ss) 

UFEM,  donor 
(hh:mm:ss) 

500 

03:46:09 

01:07:47 

00:52:45 

1000 

04:49:08 

01:28:15 

01 :04:43 

1500 

04:09:38 

01:15:55 

00:50:25 

2000 

04:23:05 

01:27:55 

01 :07:45 

2500 

03:56:44 

01:23:02 

01:03:35 

3000 

03:39:38 

01:21:03 

01:00:56 
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Figure  18:  Simulation  results  using  X-FEM  algorithm  with  projectile  velocity  3000  m/s. 
Contours  of  the  axial  velocity  are  plotted  on  the  left,  and  contours  of  the  material  number 
on  the  right.  Simulation  times  are  12,  30,  50  and  100  ps. 

Figure  18  presents  a  series  of  views  of  the  progression  of  the  formation  of  the  penetration 
channel  and  stagnation/recoil  of  the  penetrator,  using  ALEGRA  with  the  X-FEM  algorithm  for 
the  projectile  velocity  of  3000  m/s.  Four  times  are  presented  with  the  left  side  of  each  panel 
plotting  the  z-velocity  field,  and  the  right  side  showing  the  local  material  number.  The  input 
decks  for  the  simulations  are  provided  in  Appendix  1.  We  observe  that  using  X-FEM,  the 
penetrator  detaches  from  the  end  of  the  channel  after  stagnation.  This  recoil  effect  was  not  seen 
in  any  U-FEM  simulations,  which  force  the  materials  to  remain  joined. 
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4.4.  Model  boundary  condition  and  geometry  effect 


To  examine  the  model  boundary  condition  and  geometry  effect,  the  test  cases  described  in  Table 
1  and  Table  6  are  used  for  the  X-FEM  simulation.  Figure  19  and  Figure  20  show  the  meshes  for 
Test  cases  4xl8fixed  and  5xl2fixed,  respectively.  The  cell  sizes  of  the  meshes  are  constructed  as 
close  as  possible  to  the  baseline  cell  size  to  reduce  the  difference  due  to  the  cell  size.  The  input 
deck  for  Test  cases  4xl8fixed  and  5xl2fixed  are  provided  in  Appendices  2  and  3,  respectively. 


Table  6:  Comparison  of  mesh  and  boundary  for  each  test  case. 


Test  case 

Type 

Input 

cell 

size 

(mm) 

Cell  size  in 
penetration  area 
(mm) 

No.  of 
elements 

Boundary  Condition 

Figure 

Ar 

Az 

h 

Center 

Edge 

Back 

Baseline 

Graded 

0.15 

0.052 

0.073 

0.063 

301644 

Fixed 

Fixed 

Fixed 

Figure  10 

4x12free 

Graded 

0.15 

0.052 

0.073 

0.063 

301644 

Fixed 

Fixed 

Free 

Figure  10 

4x18fixed 

Graded 

0.15 

0.059 

0.075 

0.067 

317955 

Fixed 

Fixed 

Fixed 

Figure  19 

5x12fixed 

Graded 

0.15 

0.048 

0.083 

0.065 

294000 

Fixed 

Fixed 

Fixed 

Figure  20 
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Figure  19:  Computational  mesh  used  to  examine  target  plate  thickness  (Unit  on  axes:  m). 
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Figure  20:  Computational  mesh  used  to  examine  the  target  plate  radius.  (Unit  on  axes:  m). 


4.4.1.  Target  boundary 

To  represent  the  semi-infinite  impact/penetration  scenario,  no  displacement  boundary  condition 
is  applied  on  the  center,  right,  and  back  of  the  plate  in  the  baseline  model.  However,  the  steel 
target  cannot  be  semi-infinite  in  the  actual  experimental  environment,  i.e.  the  thickness  of  target 
plate  cannot  be  infinite.  To  represent  the  experimental  environment,  the  free  boundary  condition 
can  be  applied  on  the  plate  back.  Which  boundary  condition  is  better  expression  for  the  actual 
environment  is  studied  in  this  section. 

The  normalized  penetrations  from  the  models  with  the  plate  back  fixed  (Baseline  model)  and 
plate  back  free  boundary  conditions  are  slightly  different  as  shown  Figure  21  and  Table  7.  Figure 
22  shows  the  predicted  DOP,  rod  head,  and  plate  surface  histories  from  the  model  with  plate 
back  free  boundary  condition.  The  rod  head  depth  increases  continually  with  time,  after  a  brief 
drop  just  past  the  peak  as  shown  Figure  15  in  the  baseline  case.  We  speculate  that  the  motion  of 
the  plate  after  projectile  stagnation  in  the  free -boundary  case  increases  the  DOP  in  a  way  that 
would  not  occur  in  a  truly  semi-infinite  target.  To  represent  the  experimental  target  plate,  which 
was  effectively  semi-infinite,  it  seems  that  the  fixed  boundary  condition  is  appropriate.  For  the 
projectile  velocity  of  500  m/s,  the  X-FEM  simulation  stopped  at  19.9  ps,  reporting  an 
“impossible  node  position”  in  the  remap  algorithm. 
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Figure  21 :  Comparison  of  the  normalized  penetration  depths  from  the  models  with  plate 
back  fixed  (Baseline)  and  plate  back  free  boundary  conditions. 


Table  7:  Predicted  OOP  and  P/L  from  the  model  with  plate  back  fixed  and  free  boundary 
conditions  with  polynomial-fit  experimental  data.  (The  bold  fonts  indicate  the  closer 
result  to  the  experimental  data). 


Projectile 

Velocity 

(m/s) 

X-FEM,  Back 
Fixed 

X-FEM,  Back 
Free 

Data  (fit) 

DOP 

(mm) 

P/L 

DOP 

(mm) 

P/L 

DOP 

(mm) 

P/L 

500 

1.98 

0.040 

1.97 

0.039 

2.07 

0.041 

1000 

17.18 

0.344 

16.83 

0.337 

18.55 

0.371 

1500 

44.55 

0.891 

44.14 

0.883 

45.42 

0.908 

2000 

63.47 

1.269 

63.54 

1.271 

63.82 

1.276 

2500 

71.95 

1.439 

72.09 

1.442 

72.06 

1.441 

3000 

76.31 

1.526 

77.18 

1.544 

73.75 

1.475 
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Figure  22:  Predicted  DOP,  rod  head  and  plate  surface  histories  with  six  velocities  of  rod 
from  the  model  with  plate  back  free  boundary  condition. 


4.4.2.  Target  plate  thickness 

To  investigate  how  the  target  plate  thickness  in  the  model  affects  the  result,  an  X-FEM 
simulation  is  conducted  with  the  mesh  as  shown  in  Figure  19.  The  target  plate  thickness  is 
extended  to  180  mm,  and  its  radius  is  kept  at  40  mm.  The  other  conditions  such  as  material 
properties  and  boundary  conditions  are  the  same  as  the  baseline.  The  normalized  penetration 
depths  for  six  velocities  from  Test  Cases  4xl2fixed  (baseline)  and  4xl8fixed  are  close  to  each 
other  as  shown  Figure  23.  The  predicted  DOP  and  P/L  from  the  baseline  model  are  slightly 
closer  to  the  experimental  data  as  listed  in  Table  8.  In  conclusion,  the  target  plate  thickness  of 
120  mm  is  enough  to  represent  the  experimental  environment. 

Figure  24  shows  the  predicted  DOP,  rod  head,  and  plate  surface  histories  from  the  model  with 
the  target  plate  thickness  of  180  mm.  The  results  are  almost  the  same  as  the  baseline’s  (Figure  15) 
except  the  depth  of  rod  head  is  slightly  deeper  for  the  high  velocity  projectile  (3000  and  2500 
m/s).  The  X-FEM  simulation  stopped  at  16.9  ps  and  62.5  ps  when  the  projectile  velocities  are 
500  m/s  and  1000  m/s,  respectively,  again  reporting  an  “impossible  node  position”  in  the  remap 
operation. 
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1.6 


Figure  23:  Comparison  of  the  normalized  penetration  depths  from  the  models  with  plate 
thicknesses  of  12  cm  (Baseline)  and  18  cm. 

Table  8:  Predicted  OOP  and  P/L  from  the  modei  with  target  plate  thicknesses  of  120  mm 
and  180  mm  with  poiynomiai-fit  experimental  data.  (The  bold  fonts  indicate  the  closer 
result  to  the  experimental  data). 


Projectile 

Velocity 

(m/s) 

4x12  (Baseline) 

4x18fixed 

Data  (fit) 

DOP 

(mm) 

P/L 

DOP 

(mm) 

P/L 

DOP 

(mm) 

P/L 

500 

1.98 

0.040 

1.96 

0.039 

2.07 

0.041 

1000 

17.18 

0.344 

17.07 

0.341 

18.55 

0.371 

1500 

44.55 

0.891 

44.40 

0.888 

45.42 

0.908 

2000 

63.47 

1.269 

63.39 

1.268 

63.82 

1.276 

2500 

71.95 

1.439 

72.28 

1.446 

72.06 

1.441 

3000 

76.31 

1.526 

77.21 

1.544 

73.75 

1.475 
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Figure  24:  Predicted  DOP,  rod  head  and  plate  surface  histories  with  six  veiocities  of  rod 
from  the  model  with  target  plate  thickness  of  180  mm. 


4.4.3.  Target  plate  radius 

To  investigate  how  the  target  plate  radius  in  the  model  affects  the  result,  an  X-FEM  simulation  is 
conducted  with  the  mesh  as  shown  in  Figure  20  (Target  plate  thickness  is  120  mm,  radius  is  50 
mm).  The  other  conditions  such  as  material  properties,  boundary  conditions,  etc.  are  the  same  as 
the  baseline.  The  normalized  penetration  depths  for  six  velocities  from  Test  Case  4xl2fixed 
(baseline)  and  5xl2fixed  are  close  each  other  as  shown  Figure  25  and  Table  9.  In  conclusion,  the 
target  plate  radius  of  40  mm  is  enough  to  represent  the  experimental  environment. 

Figure  26  shows  the  predicted  DOP,  rod  head,  and  plate  surface  histories  from  the  model  with 
the  target  plate  radius  of  50  mm.  The  results  are  almost  the  same  as  the  baseline’s  (Figure  15). 
The  X-FEM  simulation  stopped  at  27.7  ps  and  65.6  ps  when  the  projectile  velocities  are  500  m/s 
and  1000  m/s,  respectively.  The  error  messages  reported  “impossible  node  position”  in  remap  for 
the  500-m/s  case,  and  “significant  overfill”  in  the  1000-m/s  case. 
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Figure  25:  Comparison  of  the  normalized  penetration  depths  from  the  models  with  plate 
radius  of  40  mm  (Baseline)  and  50  mm. 


Tabie  9:  Predicted  OOP  and  P/L  from  the  modei  with  target  plate  radius  of  40  mm  and  50 
mm  with  experimental  data.  (The  bold  fonts  indicate  the  cioser  result  to  the  experimental 
data).  _ _ _ _ 


Projectile 

Velocity 

(m/s) 

4x12fixed 

(Baseline) 

5x12fixed 

Experimental 

Data 

DOP 

(mm) 

P/L 

DOP 

(mm) 

P/L 

DOP 

(mm) 

P/L 

500 

1.98 

0.040 

1.99 

0.040 

2.05 

0.041 

1000 

17.18 

0.344 

17.20 

0.344 

19.68 

0.393 

1500 

44.55 

0.891 

44.27 

0.885 

45.67 

0.912 

2000 

63.47 

1.269 

62.93 

1.259 

64.75 

1.293 

2500 

71 .95 

1.439 

71.83 

1.437 

72.91 

1.456 

3000 

76.31 

1.526 

76.27 

1.525 

74.21 

1.482 
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Figure  26:  Predicted  DOP,  rod  head  and  plate  surface  histories  with  six  veiocities  of  rod 
from  the  model  with  target  plate  radius  of  50  mm. 


5.  SUMMARY  AND  CONCLUSIONS 


The  simulations  conducted  in  this  study  provide  the  first  evaluation  of  the  X-FEM 
implementation  in  ALEGRA  for  classical  problems  in  high-velocity  impact  and  penetration.  The 
results  demonstrate  that  the  new  algorithm  reproduces  the  excellent  results  obtained  using  the 
standard  algorithms  in  AEEGRA,  and  retains  their  accuracy.  Depths  of  penetration  computed 
using  X-FEM  over  the  impact  velocity  range  500  m/s  <  Vp  <  3000  m/s  differ  from  the 
experimental  data  of  Hohler  and  Stilp  (1977,1987)  [8,9]  by  only  1.6%.  They  differ  from  the 
standard-formulation  results  with  high-order  remap  by  only  0.5%.  They  also  provide  improved 
depth  of  penetration  compared  to  the  standard-formulation  results  with  low-order  remap  (see 
Table  4  and  Eigures  17  and  18). 

The  various  sensitivities  of  the  computed  X-FEM  solutions  have  been  examined  here.  The 
solution  converges  at  a  rate  that  is  consistent  with  the  expected  behavior  for  solutions  containing 
shocks:  the  depth  of  penetration  converges  smoothly  and  monotonically  at  first  order 
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(convergence  rate  1.057  <  c  <  1.066,  see  Figure  14).  Sensitivities  to  domain  size  and  boundary 
conditions  for  this  particular  problem  have  also  been  examined  to  ensure  correct  setup  of  the 
problem.  Sensitivity  to  these  parameters  is  relatively  weak,  so  that  the  baseline  configuration 
used  here  is  shown  to  be  adequate  (see  Tables  7-9  and  associated  figures). 

The  results  here  show  that  the  performance  cost  of  the  upgrade  from  standard  methods  to  X- 
FEM  is  a  factor  of  3  to  4  in  compute  time  for  this  class  of  problems  (see  Table  5).  The 
implementation  has  been  tested  for  parallel  scaling  with  good  results.  However,  no  concerted 
attempt  has  been  made  to  optimize  the  X-FEM  code  for  speed  and  efficient  memory  use.  This  is 
slated  for  work  in  the  near  future,  and  we  therefore  expect  improved  performance  to  be  observed 
soon.  We  also  note  that  out  of  29  total  X-EEM  simulations  conducted  in  the  present  study,  5 
stopped  before  completion  with  errors,  as  described  in  Section  4.  The  robustness  of  the 
implementation  is  expected  to  improve  as  current  active  development  efforts  proceed.  A  3D 
extension  of  the  X-EEM  capability  is  also  in  progress,  and  we  anticipate  testing  of  the  3D 
capability  using  the  Hohler-Stilp  data  in  future  work  as  well. 

The  ballistic  environment  considered  here  is  dominated  by  hydrodynamic  interaction  of  the 
penetrator  and  target  materials.  The  present  results  show  that  the  X-EEM  reproduces  the 
accuracy  of  the  standard  Eulerian  method  in  AEEGRA  that  has  already  been  established  for  this 
class  of  problems.  The  high  quality  of  the  results  produced  using  AEEGRA  with  X-EEM  for 
these  cases  suggests  that  this  software  is  ready  for  production  use,  and  for  testing  in  other 
situations  where  the  accuracy  of  the  standard  Eulerian  method  has  not  been  well  established. 
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APPENDIX:  INPUT  DECK 


[units.txt] 

$Um't  conversion  (A_B  =  A  to  B)  : 

$ 

$Angular: 

$  rad  =  {rad_deg  =180/3.1415926536}  deg 
$  deg  =  {deg_rad  =l/rad_deg}  rad 
$ 

Stength : 

$  mm  =  {mm_m=0 . 001  }  m 

$  m  =  {m_mm=l/mm_m}  mm 

$  cm  =  {cm_mm=10  }  mm 

$  mm  =  {mm_cm=l/cm_mm}  cm 
$  cm  =  {cm_m=0.01  }  m 

$  m  =  {m_cm=l/cm_m}  cm 

$  ft  =  {ft_m=0.3048}  m 

$  m  =  {m_ft=l/ft_m}  ft 

$  in  =  {in_m=0.0254}  m 

$  m  =  {m_i n=l/i n_m}  i n 

$ 

$Pressure : 

$  MPa  =  {MPa_Pa  =  1E6}  Pa 

$  Pa  =  {Pa_MPa  =  l/MPa_Pa}  MPa 

$  psi  =  {psi_Pa=6894.757}  Pa 

$  Pa  =  {Pa_psi=l/psi_Pa}  psi 

$  Pa  =  {Pa_psf=0. 0208854}  psf 

$  psf  =  {psf_Pa=l/Pa_psf}  Pa 

$  MPa  =  {MPa_psf=MPa_Pa*Pa_psf}  psf 

$  MPa  =  {MPa_psi=MPa_Pa*Pa_psi}  psi 

$  psf  =  {psf_psi=psf_Pa*Pa_psi}  psi 

$ 

$Ti me : 


$ 

mus 

= 

{mus_s 

= 

le-6 

} 

s 

$ 

mi  n 

= 

{mi n_s 

= 

60 

} 

s 

$ 

h 

= 

{h  min 

= 

60 

} 

mi  n 

$ 

d 

= 

{d_h 

= 

24 

} 

h 

$ 

mon 

{mon_d 

30.4166666667} 

d 

$ 

yr 

= 

{yr_d 

= 

365 

} 

d 

$ 

dec 

= 

{dec_yr 

= 

10 

} 

yr 

$ 

cen 

= 

{cen_dec 

= 

10 

} 

dec 

$ 

mi  1 

= 

{mil  cen 

= 

10 

} 

cen 

$ 

h 

= 

{h_s 

= 

h  min*min  s 

} 

s 

$ 

d 

{d_s 

d_h*h_s 

} 

s 

$ 

mon 

= 

{mon_s 

= 

mon_d*d_s 

} 

s 

$ 

yr 

{yr_s 

yr_d*d_s 

} 

s 

$ 

dec 

= 

{dec_s 

= 

dec_yr*yr_s 

} 

s 

$ 

cen 

= 

{cen_s 

= 

cen_dec*dec_ 

-S} 

s 

$ 

mi  1 

= 

{mi  l_s 

= 

mi  l_cen*cen_ 

-S} 

s 

$ 

s 

{s_mus 

l/mus_s 

} 

mus 

$ 

s 

= 

{s  min 

= 

1/min  s 

} 

mi  n 

$ 

s 

= 

{s_h 

= 

l/h_s 

} 

h 

$ 

s 

= 

{s_mon 

= 

l/mon_s 

} 

mon 

$ 

s 

{s_yr 

l/yr_s 

} 

yr 

Appendix  1 :  Input  deck  for  4x12fixed 


[parameters.txt] 

$ 

$  Define  Parameters 
$ 


$  Test  Case:  T_CAS  =  {T_CAS='subM'} 

$  Cell  size:  S_CEL  =  {S_CEL=0.15*mm_m}  m 

$  Model  Center  (X):  C_MDL  =  {C_MDL=  0.00}  m 

$  Radius  of  Plate  (X):  R_PLT  =  {R_PLT=  40.0*mm_m}  m 

$  Thickness  of  Plate  (Y) :  T_PLT  =  {T_PLT=  120.0*mm_m}  m 

$  Radius  of  Rod  (X):  R_ROD  =  {R_ROD=  2.5*mm_m}  m 

$  Length  of  Rod  (Y) :  L_ROD  =  {L_ROD=  50.0*mm_m}  m 

$  Velocity  of  Rod:  V_ROD  =  {V_ROD=  500}  m/s 

$  Front  of  Plate  (Y) :  F_PLT  =  {F_PLT=  0.00}  m 

$  Back  of  Plate  (Y) :  B_PLT  =  {B_PLT=  F_PLT+T_PLT}  m 

$  Head  of  Rod  (Y) :  H_ROD  =  {H_ROD=  F_PLT}  m 

$  Tail  of  Rod  (Y) :  T_ROD  =  {T_ROD=  H_ROD-L_ROD}  m 

$  Back  of  Model  (Y) :  B_MDL  =  {B_MDL=  R_PLT}  m 

$  Front  of  1st  Graded  Area  (Y) :  F_1GA  =  {F_lGA=  78.*mm_m}  m  (hole  depth  is  0.0741  at  3000  m/s) 
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$  Radius  of  1st  Graded  Area  (X):  R_lGA  =  {R_lGA=  4.*R_ROD}  m 
$  Gap  between  Rod  tail  and  Front  of  Model  (Y) :  L_GAP  =  {L_GAP=  S_CEL}  m 
$  Front  of  Model  (Y) :  F_MDL  =  {F_MDL=  T_ROD-L_GAP}  m 

$  Length  of  Model  (Y) :  L_MDL  =  {L_MDL=  T_PLT+L_ROD+L_GAP}  m 

Note:  V_ROD  varied  from  500  m/s  to  3000  m/s  for  the  projectile  rod  velocities. 


[hohlerstilprha.inp]  for  subcase  5-G 


$  {i ncl ude("/home/bypark/common/uni ts . txt")} 

$  {_FORMAT="%.10g"} 

$  -  COMMENTS  - 

$  Author:  Brian  Leavy 

$  Modified:  Byoung  Yoon  Park  12-02-2011 

$  Experiment:  Hohler  &  Stilp  W  rod  vs  RHA 

$  Comments:  Validation  test  for  X-FEM  in  ALEGRA 

$===  VARIABLES==== 

$  {includeC'parameters.txt")} 

$  Time  Interval:  {T_INT  =  l.Oe-6}  s 

$  Plot  Interval:  {P_INT  =  T_INT}  s 

$  ===  execution  control  === 

$ - Initiation  and  Termination 

Title 

Hohler-Stilprha  {T_CAS},  v={V_ROD}  m/s,  R={R_ROD*m_mm}  mm,  L={L_ROD*m_cm}  cm,  CELL={S_CEL*m_mm} 
mm 

units,  si 

READ  RESTART  DUMP  =  -1  $  Restart  at  the  latest  available  restart  dump 

termination  time,  {tt=160.0e-6}  $  Simulation  terminated  at  {tt*s_mus}  mus. 

TERMINATION  CPU  {tc=6.*h_s}  $  Run  will  be  terminated  within  CPU  time  lhr/64nodes(Agg) 

Stermi nation  cycle,  1 

CHECK  SHUTDOWN  FILE  {csf=2}m  $  Check  shutdown  file  every  {csf}  minutes 

EMIT  RESTART,  TIME  INTERVAL  {T_INT}  $  Make  dmp  file  every  {T_lNT*s_mus}  mus. 

$ -  I/O  Control 

emit  plot,  time  interval  ={P_INT}  $  s 

$emit  plot,  cycle  interval=10  from  time  7.e-5  to  15.e-5 

emit  hisplt,  cycle  interval=10 

$emit  hisplt,  time  i nterval={T_lNT}  $  s 

$emit  hisplt,  cycle  interval=l  from  time  7.e-5  to  15.e-5 

plot  variables 
no  underscores 
velocity 
density 
density, avg 
temperature 
temperature, avg 
pressure 
pressure , avg 
mat_mi n_coords 
mat_max_coo  rd  s 
end 

history  plot  variables 
no  underscores 
velocity 
density 
temperature 
pressure 
end 


spy 

PlotTime(0.0,{P_lNT})  ; 

define  main() 

\{ 

pprintfC"  PLOT:  Cycle=%d,  Time=%e\n" , CYCLE, TIME) ; 

%  ListVariablesC)  ; 

Image ("Hohl ersti 1 pRha" , WHITE , BLACK) ; 

Label ("{T_CAS} ,  v={V_ROD}  m/s.  Graded  cell={S_CEL*m_mm}  mm"); 
SMOOTH_SHADING  =  OFF; 

ColorMapRangeCO. ,20000. ,LIN_MAP) ; 

%  AutoScaleMap("DENSlTY",0. ,20000.) ; 

Drawcol orMap ("DENSITY" ,0.1, 0.3, 0.3, 0.8); 
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Plot2D("DENSITY") ; 

Color(PURPLE)  ; 

Draw2DMatContour ; 

Color(BLACK) ; 

DrawText(spn' ntf ("t  =  %.2f  ~m~s",l.e6*TlME)  ,0.4,0.05)  ; 

Color(BLACK) ; 

Logoshadi ng(ON) ; 

LogoText(ON) ; 

LogoColors(BLACK, CYAN, BLACK) ; 

LogoPosi ti on (0 . 0 , 0 . 0 , 0 . 4 , 0 . 15) ; 

DrawSNL_Logo(l. 3) ; 

Endimage ; 

\} 

endspy 

$  ===  PHYSICS  INPUT=== 
solid  dynamics 
cyl  i  ndri  cal 

constant  volume  fraction  algorithm 
end 

xfem 

tangential  search  tolerance,  0.05 
maxenf  iterations,  100 
gap  rate  tolerance,  l.e-8 
mspair  2,1 
end 

domai n 

remap  method,  intersection 
pattern  reconstruction 
interface  order  201 
pattern  smoothing  =  off 
pattern  fragment  motion  =  off 
end 

time  step  scale,  {tss=.4}  $  {tss}  times  of  the  value  of  the  calculated  stable  time  step. 

$ - MESH - 

$BX  :{BX  =  1}  $  Number  of  Blocks  in  X 
$BY  :{BY  =  1}  $  Number  of  Blocks  in  Y 
$NX_TOT  :{NX_TOT  =  i nt ( (r_plt/s_CEL)/bx) } 

$NY_TOT  :{NY_TOT  =  i nt ( (L_MDL/S_CEL)/BY) } 

$NX  : {NX  =  NX_TOT/BX} 

$NY  : {NY  =  NY_TOT/BY} 

SNSLOTS  :{NSL0TS  =  BX*BY} 

$NUM_CELLS  :{NUM_CELLS  =  NX_TOT*NY_TOT} 

$R_DER  :{R_DER  =  4.*R_ROD}  $  Dense  element  radius  bound  (R-di r) 

$Z_DEU  :{Z_DEU  =  0.078}  $  Dense  element  upper  bound  (Z-di r)  (hole  depth  is  0.0741  at  3000  m/s) 

$Z_DEL  :{Z_DEL  =  0.000}  $  Dense  element  lower  bound  (Z-di r) 

$F_DER  :{F_DER  =  0.1}  $  Increase  Rate  (R-di r) ,  The  larger,  the  smaller  interval 

$F_DER  :{F_DEU  =  0.7}  $  Increase  Factor  (+Z-di r) ,  The  larger,  the  smaller  interval 

$F_DEL  :{F_DEL  =  0.3}  $  Increase  Factor  (-Z-di r) ,  The  larger,  the  larger  interval 

mesh,  inline 
recti  1 i near 
bx  =  {BX} 
by  =  {BY} 
nx  =  {NX} 
ny  =  {NY} 

gmin  =  {C_MDL}  {F_MDL} 

gmax  =  {NX_TOT*BX*S_CEL}  {F_MDL+NY_TOT*BY*S_CEL} 
end 

set  assign 

sideset,ilo, 30 
si  deset, i hi ,10 
sideset, qlo,40 
sideset,]hi ,20 
end 

USER  DEFINED  ELEMENT  DENSITY,  X 
fi  el  d=l; 

if(coord  <  {R_DER})  {"{"}field  =  1;{"}"} 

if(coord  >=  {R_DER})  {"{"}field  =  l.*exp(-100.*((coord)-{R_DER*F_DER})) ;{"}"} 

END 

USER  DEFINED  ELEMENT  DENSITY,  Y 
fi  el  d=l; 
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if(coord  <  {Z_DEU})  {"{"}field  = 

ifCcoord  >=  {Z_DEU})  {"{"}field  =  l.*exp(-80.*((coord)-{Z_DEU*F_DEU})) ;{ 
if(coord  <=  {Z_DEL})  {"{"}field  =  l.*exp(  70.*((coord)-{Z_DEU*F_DEL})) ;{ 
ifCcoord  <  {T_ROD})  {"{"}field  =  !;{"}''} 


END 

end 

$$ 

$  -  BOUNDARY  CONDITION  - 

no  displacement,  sideset  20, z  $  +Z 
no  displacement,  sideset  10, r  $  +R 
no  displacement,  sideset  30, r  $  -R  =  0 

block  1 

EULERIAN  MESH 
add  diatom  input 
end 


$  -  MATERIAL  INSERTION  - 

$ - 

$  material  1,  rod,  93%  WHA,  Tungsten 
$  material  2,  plate,  HzB  Steel 

$ - 

di atom 

translateCO. ,{H_ROD}) 
package  'ROD' 
material  1 
yvelocity  {V_ROD} 
insert  box 

pi  {C_MDL}  {H_ROD} 

p2  {R_ROD}  {T_ROD} 

endi 
endp 
endt 

translateCO. ,{F_PLT}) 
package  'PLATE' 
material  2 
insert  box 

pi  {C_MDL}  {F_PLT} 

p2  {R_PLT}  {B_PLT} 

endi 
endp 
endt 

enddiatom 


$  -  TRACERS  - 

tracer  points 

lagrangian  Tracer  1  x 
lagrangian  tracer  2  x 
lagrangian  tracer  3  x 
lagrangian  tracer  4  x 
lagrangian  tracer  5  x 
end 

plot,  exotracer 

file  =  'tracer.exo' 
end 

$  plot,  vtk  mesh 
$  file  =  'hohlerstilprha 
$  end 


{C_MDL}  y  {H_ROD} 
{C_MDL}  y  {T_ROD} 
{C_MDL}  y  {F_PLT} 
{R_PLT-S_CEL}  y  {F_PLT} 
{R_ROD}  y  {F_PLT} 


$  Center  of  Rod  Head 

$  Center  of  Rod  Tail 

$  Center  of  Plate  Front 

$  Edge  of  Plate  Front  - 

$  Edge  of  Rod  on  Plate 


end 


$  ===  MATERIAL  === 
material  1  $  93%  WHA 

model  =  1 
end 

model  1  cth  elastic  plastic 
eos  model  =  100 

yield  model  =  10 

$  fracture  model  =  1000 
end 

model  10  johnson  cook  ep 

ajo  =  1.365e9  $  Hohler  &  Stilp 

bjo  =  0.1765e9 

cjo  =  0.016 

mjo  =  1.00 

njo  =  0.12 

tjo  =  3695. 

poisson  =  0.281 


-  Tracer  H 

-  Tracer  P 
Front 


43 


$  TUNGSTEN 

$  Scaled  for  17.760  g/cc 


end 

model  100,  keos  sesame 
feos  =  'sesame' 
neos  =  3550 
sr  =  1.094 
clip  =  1.0 
end 

material  2  $  HzB  Steel 

model  =  2 
end 

model  2  cth  elastic  plastic 
eos  model  =  200 

yield  model  =  20 

end 

model  20  johnson  cook  ep 

ajo=0.810e9  $  Hohler  &  Stilp  HzB 

bjo=0.5095e9 
cj 0=0. 014 
mj 0=1. 030 
njo=0 . 260 
tj 0=1818 . 
poisson  =  0.299 
end 

model  200,  keos  sesame 
feos  =  'sesame' 
neos  =  2150  $  Iron 

sr  =  1.003  $  Scaled  for  7.85  g/cc 

clip  =  1.0 
end 

exit 


[hohlerstilprha.inp]  for  U-FEM  aggressive  remap  simulation 


$  {i ncl ude("/home/bypark/common/uni ts . txt")} 

$  {_FORMAT="%.10g"} 

$  -  COMMENTS  - 

$  Author:  Brian  Leavy 

$  Modified:  Byoung  Yoon  Park  12-02-2011 

$  Experiment:  Hohler  &  Stilp  W  rod  vs  RHA 

$  Comments:  Validation  test  for  X-FEM  in  ALEGRA 

$===  VARIABLES==== 

$  {includeC'parameters.txt")} 

$  Time  Interval:  {T_INT  =  l.Oe-6}  s 

$  Plot  Interval:  {P_INT  =  T_INT}  s 

$  ===  EXECUTION  CONTROL  === 

$ - Initiation  and  Termination 

Title 

Hohler-Stilprha  {T_CAS},  v={V_ROD}  m/s,  R={R_ROD*m_mm}  mm,  L={L_ROD*m_cm}  cm,  CELL={S_CEL*m_mm} 
mm 

units,  si 

READ  RESTART  DUMP  =  -1  $  Restart  at  the  latest  available  restart  dump 

termination  time,  {tt=160 . Oe-6}  $  Simulation  terminated  at  {tt*s_mus}  mus. 

TERMINATION  CPU  {tc=6.*h_s}  $  Run  will  be  terminated  within  CPU  time  lhr/64nodes(Agg) 

Stermi nation  cycle,  1 

CHECK  SHUTDOWN  FILE  {csf=2}m  $  Check  shutdown  file  every  {csf}  minutes 

EMIT  RESTART,  TIME  INTERVAL  {T_INT}  $  Make  dmp  file  every  {T_lNT*s_mus}  mus. 

$ -  I/O  Control 

emit  plot,  time  interval  ={P_INT}  $  s 
emit  hisplt,  cycle  interval=10 

plot  variables 
no  underscores 
velocity 
density 
density, avg 
temperature 
temperature , avg 
pressure 
pressure , avg 
mat_mi n_coords 
mat_max_coo  rd  s 
end 
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history  plot  variables 
no  underscores 
velocity 
density 
temperature 
pressure 
end 


spy 

PlotTime(0.0,{P_lNT})  ; 

define  main() 

\{ 

pprintfC"  PLOT:  Cycle=%d,  Time=%e\n" .CYCLE, TIME) ; 

%  ListVariablesC)  ; 

Image ("Hohl ersti 1 pRha" .WHITE , BLACK) ; 

Label ("{T_CAS} ,  v={V_ROD}  m/s,  Graded  cell={S_CEL*m_mm}  mm"); 

SMOOTH_SHADING  =  OFF; 

ColorMapRange(0. ,20000. ,LIN_MAP) ; 

%  AutoScaleMapC'DENSiTY"  ,0.  ,20000.  )  ; 

Drawcol orMap ("DENSITY" ,0.1, 0.3, 0.3, 0.8); 

Plot2D("DENSITY") ; 

Color(PURPLE)  ; 

Draw2DMatContour ; 

Color(BLACK) ; 

DrawText(spri ntf ("t  =  %.2f  ~m~s",l.e6*TlME) ,0.4,0.05) ; 

Color(BLACK) ; 

Logoshadi ng(ON) ; 

LogoText(ON) ; 

LogoColors(BLACK, CYAN, BLACK) ; 

LogoPosi ti on (0 . 0 , 0 . 0 , 0 . 4 , 0 . 15) ; 

DrawSNL_Logo(l. 3) ; 

Endimage ; 

\} 

endspy 

$  ===  PHYSICS  INPUT=== 
solid  dynamics 
cyl  i  ndri  cal 

constant  volume  fraction  algorithm 
end 

domai n 

pattern  reconstruction 
interface  order  201 
pattern  smoothing  =  off 
pattern  fragment  motion  =  off 
end 

time  step  scale,  {tss=.4}  $  {tss}  times  of  the  value  of  the  calculated  stable  time  step. 

$ - MESH - 

$BX  :{BX  =  1}  $  Number  of  Blocks  in  X 
$BY  :{BY  =  1}  $  Number  of  Blocks  in  Y 
$NX_TOT  :{NX_TOT  =  i nt ( (r_plt/s_CEL)/bx) } 

$NY_TOT  :{NY_TOT  =  i nt ( (L_MDL/S_CEL)/BY) } 

$NX  : {NX  =  NX_TOT/BX} 

$NY  : {NY  =  NY_TOT/BY} 

SNSLOTS  :{NSL0TS  =  BX*BY} 

$NUM_CELLS  :{NUM_CELLS  =  NX_TOT*NY_TOT} 

$R_DER  :{R_DER  =  4.*R_ROD}  $  Dense  element  radius  bound  (R-di r) 

$Z_DEU  :{Z_DEU  =  0.078}  $  Dense  element  upper  bound  (Z-di r)  (hole  depth  is  0.0741  at  3000  m/s) 

$Z_DEL  :{Z_DEL  =  0.000}  $  Dense  element  lower  bound  (Z-di r) 

$F_DER  :{F_DER  =0.1}  $  Increase  Rate  (R-di r) ,  The  larger,  the  smaller  interval 

$F_DER  :{F_DEU  =  0.7}  $  Increase  Factor  (+Z-di r) ,  The  larger,  the  smaller  interval 

$F_DEL  :{F_DEL  =  0.3}  $  Increase  Factor  (-Z-di r) ,  The  larger,  the  larger  interval 

mesh,  inline 
recti  1 i near 
bx  =  {BX} 
by  =  {BY} 
nx  =  {NX} 
ny  =  {NY} 

gmin  =  {C_MDL}  {F_MDL} 

gmax  =  {NX_TOT*BX*S_CEL}  {F_MDL+NY_TOT*BY*S_CEL} 
end 
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set  assign 

sideset.ilo, 30 
si  deset, i hi ,10 
sideset, j1o,40 
sideset,]hi ,20 
end 

USER  DEFINED  ELEMENT  DENSITY,  X 
fi  el  d=l; 

if(coord  <  {R_DER})  {"{"}field  = 

if(coord  >=  {R_DER})  {"{"}field  =  l.*exp(-100.*((coord)-{R_DER*F_DER})) ;{"}"} 

END 

USER  DEFINED  ELEMENT  DENSITY,  Y 
fi  el  d=l; 

if(coord  <  {z_DEU})  {"{"}field  = 

ifCcoord  >=  {Z_DEU})  {"{''}field  =  l.*exp(-80.*((coord)-{Z_DEU*F_DEU}))  ;{"}"} 
ifCcoord  <=  {Z_DEL})  {"{"}field  =  l.*exp(  70. *((coord)-{Z_DEU*F_DEL})) ;{"}"} 
if(coord  <  {T_ROD})  {"{"}field  = 

END 

end 

$$ 

$  -  BOUNDARY  CONDITION  - 

no  displacement,  sideset  20, z  $  +Z 
no  displacement,  sideset  10, r  $  +R 
no  displacement,  sideset  30, r  $  -R  =  0 

block  1 

EULERIAN  MESH 
add  diatom  input 
AGGRESSIVE  ADVECTION 
end 

$  -  MATERIAL  INSERTION  - 

$ - 

$  material  1,  rod,  93%  WHA,  Tungsten 
$  material  2,  plate,  HzB  Steel 

$ - 

di atom 

translate(0. ,{H_ROD}) 
package  'ROD' 
material  1 
yvelocity  {V_ROD} 
insert  box 

pi  {C_MDL}  {H_ROD} 

p2  {R_ROD}  {T_ROD} 

endi 
endp 
endt 

translate(0. ,{F_PLT}) 
package  'PLATE' 
material  2 
insert  box 

pi  {C_MDL}  {F_PLT} 

p2  {R_PLT}  {B_PLT} 

endi 
endp 
endt 

enddiatom 

$  -  TRACERS  - 

tracer  points 

lagrangian  Tracer  1  x  {C_MDL}  y  {H_ROD}  $  Center  of  Rod  Head 

lagrangian  tracer  2  x  {C_MDL}  y  {T_ROD}  $  Center  of  Rod  Tail 

lagrangian  tracer  3  x  {C_MDL}  y  {F_PLT}  $  Center  of  Plate  Front 

lagrangian  tracer  4  x  {R_PLT-S_CEL}  y  {F_PLT}  $  Edge  of  Plate  Front 

lagrangian  tracer  5  x  {R_ROD}  y  {F_PLT}  $  Edge  of  Rod  on  Plate  Front 

end 

plot,  exotracer 

file  =  'tracer.exo' 
end 

$  plot,  vtk  mesh 
$  file  =  'hohlerstilprha' 

$  end 

end 
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$  ===  MATERIAL  === 
material  1  $  93%  WHA 

model  =  1 
end 

model  1  cth  elastic  plastic 
eos  model  =  100 

yield  model  =  10 

$  fracture  model  =  1000 
end 

model  10  johnson  cook  ep 

ajo  =  1.365e9  $  Hohler  &  Stilp 

bjo  =  0.1765e9 
cjo  =  0.016 
mjo  =  1.00 
njo  =  0.12 
tjo  =  3695. 
poisson  =  0.281 
end 

model  100,  keos  sesame 
feos  =  'sesame' 

neos  =  3550  $  TUNGSTEN 

sr  =  1.094  $  Scaled  for  17.760  g/cc 

clip  =  1.0 
end 

material  2  $  HzB  Steel 

model  =  2 
end 

model  2  cth  elastic  plastic 
eos  model  =  200 

yield  model  =  20 

end 

model  20  johnson  cook  ep 

ajo=0.810e9  $  Hohler  &  Stilp  HzB 

bjo=0.5095e9 
cj 0=0. 014 
mj 0=1. 030 
njo=0 . 260 
tj 0=1818 . 
poisson  =  0.299 
end 

model  200,  keos  sesame 
feos  =  'sesame' 
neos  =  2150  $  Iron 

sr  =  1.003  $  Scaled  for  7.85  g/cc 

clip  =  1.0 
end 

exi  t 


[hohlerstilprha.inp]  for  DNR_G_P15 


$  {i ncl ude("/home/bypark/common/uni ts . txt")} 

$  {_FORMAT="%.10g"} 

$  -  COMMENTS  - 

$  Author:  Brian  Leavy 

$  Modified:  Byoung  Yoon  Park  12-02-2011 

$  Experiment:  Hohler  &  Stilp  W  rod  vs  RHA 

$  Comments:  Validation  test  for  X-FEM  in  ALEGRA 

$===  VARIABLES==== 

$  {includeC'parameters.txt")} 

$  Time  Interval:  {T_INT  =  l.Oe-6}  s 

$  Plot  Interval:  {P_INT  =  T_INT}  s 

$  ===  execution  control  === 

$ - Initiation  and  Termination 

Title 

Hohler-Stilprha  {T_CAS},  v={V_ROD}  m/s,  R={R_ROD*m_mm}  mm,  L={L_ROD*m_cm}  cm,  CELL={S_CEL*m_mm} 
mm 

units,  si 

READ  RESTART  DUMP  =  -1  $  Restart  at  the  latest  available  restart  dump 

termination  time,  {tt=160.0e-6}  $  Simulation  terminated  at  {tt*s_mus}  mus. 

TERMINATION  CPU  {tc=6.*h_s}  $  Run  will  be  terminated  within  CPU  time  lhr/64nodes(Agg) 

Stermi nation  cycle,  1 

CHECK  SHUTDOWN  FILE  {csf=2}m  $  Check  shutdown  file  every  {csf}  minutes 

EMIT  RESTART,  TIME  INTERVAL  {T_INT}  $  Make  dmp  file  every  {T_lNT*s_mus}  mus. 
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$ -  I/O  Control 

emit  plot,  time  interval ={P_INT}  $  s 
emit  hisplt,  cycle  interval=10 

plot  variables 
no  underscores 
velocity 
density 
density, avg 
temperature 
temperature, avg 
pressure 
pressure , avg 
mat_mi n_coords 
mat_max_coo  rd  s 
end 

history  plot  variables 
no  underscores 
velocity 
density 
temperature 
pressure 
end 


spy 

PlotTime(0.0,{P_lNT})  ; 

define  main() 

\{ 

pprintfC"  PLOT:  Cycle=%d,  Time=%e\n" , CYCLE, TIME) ; 

%  ListVariablesC)  ; 

Image ("Hohl ersti 1 pRha" , WHITE , BLACK) ; 

Label ("{T_CAS} ,  v={V_ROD}  m/s.  Graded  cell={S_CEL*m_mm}  mm"); 
SMOOTH_SHADING  =  OFF; 

ColorMapRangeCO. ,20000. ,LIN_MAP) ; 

%  AutoScaleMapC'DENSiTY"  ,0.  ,20000.)  ; 

DrawColorMapC'DENSlTY"  ,0. 1,0. 3 ,0. 3 ,0.8)  ; 

Plot2D("DENSITY") ; 

Color(PURPLE) ; 

Draw2DMatContour ; 

Color(BLACK) ; 

DrawText(spri ntf ("t  =  %.2f  ~m~s",l.e6*TlME) ,0.4,0.05) ; 
Color(BLACK) ; 

Logoshadi ng(ON) ; 

LogoText(ON) ; 

LogoColors(BLACK, CYAN, BLACK) ; 

LogoPosi ti on (0 . 0 , 0 . 0 , 0 . 4 , 0 . 15) ; 

DrawSNL_Logo(l. 3) ; 

Endimage ; 

\} 

endspy 


$  ===  PHYSICS  INPUT=== 
solid  dynamics 
cyl  i  ndri  cal 

constant  volume  fraction  algorithm 
end 

domai n 

pattern  reconstruction 
interface  order  201 
pattern  smoothing  =  off 
pattern  fragment  motion  =  off 
end 

time  step  scale,  {tss=.4}  $  {tss}  times  of  the  value  of  the  calculated  stable  time  step. 

$ - MESH - 

$BX  :{BX  =  1}  $  Number  of  Blocks  in  X 
$BY  :{BY  =  1}  $  Number  of  Blocks  in  Y 
$NX_TOT  :{NX_TOT  =  i nt ( (r_plt/s_CEL)/bx) } 

$NY_TOT  :{NY_TOT  =  i nt ( (L_MDL/S_CEL)/BY) } 

$NX  : {NX  =  NX_TOT/BX} 

$NY  : {NY  =  NY_TOT/BY} 
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SNSLOTS  :{NSL0TS  =  BX*BY} 

$NUM_CELLS  :{NUM_CELLS  =  NX_TOT*NY_TOT} 

$R_DER  :{R_DER  =  4.*R_ROD}  $  Dense  element  radius  bound  (R-di r) 

$Z_DEU  :{Z_DEU  =  0.078}  $  Dense  element  upper  bound  (Z-di r)  (hole  depth  is  0.0741  at  3000  m/s) 

$Z_DEL  :{Z_DEL  =  0.000}  $  Dense  element  lower  bound  (Z-di r) 

$F_DER  :{F_DER  =  0.1}  $  Increase  Rate  (R-di r) ,  The  larger,  the  smaller  interval 

$F_DER  :{F_DEU  =  0.7}  $  Increase  Factor  (+Z-di r) ,  The  larger,  the  smaller  interval 

$F_DEL  :{F_DEL  =  0.3}  $  Increase  Factor  (-Z-di r) ,  The  larger,  the  larger  interval 

mesh,  inline 
recti  1 i near 
bx  =  {BX} 
by  =  {BY} 
nx  =  {NX} 
ny  =  {NY} 

gmin  =  {C_MDL}  {F_MDL} 

gmax  =  {NX_TOT*BX*S_CEL}  {F_MDL+NY_TOT*BY*S_CEL} 
end 

set  assign 

sideset,ilo, 30 
si  deset, i hi ,10 
sideset, jlo,40 
sideset,]hi ,20 
end 

USER  DEFINED  ELEMENT  DENSITY,  X 
field=l; 

if(coord  <  {R_DER})  {"{"}field  =  !;{"}"} 

if(coord  >=  {R_DER})  {"{"}field  =  l.*exp(-100.*((coord)-{R_DER*F_DER})) ;{"}"} 

END 

USER  DEFINED  ELEMENT  DENSITY,  Y 
fi  el  d=l; 

if(coord  <  {Z_DEU})  {"{"}field  =  !;{"}"} 

if(coord  >=  {Z_DEU})  {"{"}field  =  l.*exp(-80.*((coord)-{Z_DEU*F_DEU})) ;{"}"} 
if(coord  <=  {Z_DEL})  {"{"}field  =  l.*exp(  70. *((coord)-{Z_DEU*F_DEL})) ;{"}"} 
if(coord  <  {T_ROD})  {"{"}field  =  1;{"}"} 

END 

end 

$$ 

$  -  BOUNDARY  CONDITION  - 

no  displacement,  sideset  20, z  $  +Z 
no  displacement,  sideset  10, r  $  +R 
no  displacement,  sideset  30, r  $  -R  =  0 

block  1 

eulerian  mesh 

add  diatom  input 

user  defined  advection 

mixed  cells,  donor  advection 
mixed  neighbors,  donor  advection 
mixed  stencil,  donor  advection 
element,  donor  advection 
end 
end 

$  -  MATERIAL  INSERTION  - 

$ - 

$  material  1,  rod,  93%  WHA,  Tungsten 
$  material  2,  plate,  HzB  Steel 

$ - 

di atom 

translate(0. ,{H_ROD}) 
package  'ROD' 
material  1 
yvelocity  {V_ROD} 
insert  box 

pi  {C_MDL}  {H_ROD} 

p2  {R_ROD}  {T_ROD} 

endi 
endp 
endt 

translate(0. ,{F_PLT}) 
package  'PLATE' 
material  2 
insert  box 

pi  {C_MDL}  {F_PLT} 

p2  {R_PLT}  {B_PLT} 

endi 
endp 
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endt 

enddiatom 


$  -  TRACERS  - 

tracer  points 

lagrangian  Tracer  1  x  {C_MDL}  y  {H_ROD}  $  Center  of  Rod  Head  --  Tracer  h 

lagrangian  tracer  2  x  {C_MDL}  y  {T_ROD}  $  Center  of  Rod  Tail 

lagrangian  tracer  3  x  {C_MDL}  y  {F_PLT}  $  Center  of  Plate  Front 

lagrangian  tracer  4  x  {R_PLT-S_CEL}  y  {F_PLT}  $  Edge  of  Plate  Front  --  Tracer  p 
lagrangian  tracer  5  x  {R_ROD}  y  {F_PLT}  $  Edge  of  Rod  on  Plate  Front 

end 

plot,  exotracer 

file  =  'tracer.exo' 
end 

$  plot,  vtk  mesh 
$  file  =  'hohlerstilprha' 

$  end 

end 

$  ===  MATERIAL  === 
material  1  $  93%  WHA 

model  =  1 
end 

model  1  cth  elastic  plastic 
eos  model  =  100 

yield  model  =  10 

$  fracture  model  =  1000 
end 

model  10  johnson  cook  ep 

ajo  =  1.365e9  $  Hohler  &  Stilp 

bjo  =  0.1765e9 
cjo  =  0.016 
mjo  =  1.00 
njo  =  0.12 
tjo  =  3695. 
poisson  =  0.281 
end 

model  100,  keos  sesame 
feos  =  'sesame' 

neos  =  3550  $  TUNGSTEN 

sr  =  1.094  $  Scaled  for  17.760  g/cc 

clip  =  1.0 
end 

material  2  $  HzB  Steel 

model  =  2 
end 

model  2  cth  elastic  plastic 
eos  model  =  200 

yield  model  =  20 

end 

model  20  johnson  cook  ep 

ajo=0.810e9  $  Hohler  &  Stilp  HzB 

bjo=0.5095e9 
cj 0=0. 014 
mj 0=1. 030 
njo=0 . 260 
tj 0=1818 . 
poisson  =  0.299 
end 

model  200,  keos  sesame 
feos  =  'sesame' 
neos  =  2150  $  Iron 

sr  =  1.003  $  Scaled  for  7.85  g/cc 

clip  =  1.0 
end 

exi  t 
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Appendix  2:  Input  deck  for  4x18fixed 


[parameters.txt] 

$ 

$  Define  Parameters 
$ 

$  Test  case:  T_CAS  =  {T_CAS='4xl8fx0500'} 

$  Cell  size:  S_CEL  =  {S_CEL=0.17*mm_m}  m 

$  Model  Center  (X):  C_MDL  =  {C_MDL=  0.00}  m 

$  Radius  of  Plate  (X):  R_PLT  =  {R_PLT=  40.0*mm_m}  m 

$  Thickness  of  Plate  (Y) :  T_PLT  =  {T_PLT=  180.0*mm_m}  m 

$  Radius  of  Rod  (X):  R_ROD  =  {R_ROD=  2.5*mm_m}  m 

$  Lengh  of  Rod  (Y) :  L_ROD  =  {L_ROD=  50.0*mm_m}  m 

$  Velocity  of  Rod:  V_ROD  =  {V_ROD=  500}  m/s 

$  Front  of  Plate  (Y) :  F_PLT  =  {F_PLT=  0.00}  m 

$  Back  of  Plate  (Y) :  B_PLT  =  {B_PLT=  F_PLT+T_PLT}  m 

$  Head  of  Rod  (Y) :  H_ROD  =  {H_ROD=  F_PLT}  m 

$  Tail  of  Rod  (Y) :  T_ROD  =  {T_ROD=  H_ROD-L_ROD}  m 

$  Back  of  Model  (Y) :  B_MDL  =  {B_MDL=  R_PLT}  m 

$  Front  of  1st  Graded  Area  (Y) :  F_1GA  =  {F_lGA=  78.*mm_m}  m  (hole  depth  is  0.0741  at  3000  m/s) 

$  Radius  of  1st  Graded  Area  (X):  R_1GA  =  {R_lGA=  4.*R_ROD}  m 

$  Gap  between  Rod  tail  and  Fornt  of  Model  (Y) :  L_GAP  =  {L_GAP=  S_CEL}  m 
$  Front  of  Model  (Y) :  F_MDL  =  {F_MDL=  T_ROD-L_GAP}  m 

$  Length  of  Model  (Y) :  L_MDL  =  {L_MDL=  T_PLT+L_ROD+L_GAP}  m 

Note:  V_ROD  varied  from  500  m/s  to  3000  m/s  for  the  projectile  rod  velocities. 


[hohlerstilprha.inp] 

$  {i ncl ude("/home/bypark/common/uni ts . txt") } 

$  {_FORMAT="%.10g"} 

$  -  COMMENTS  - 

$  Author:  Brian  Leavy 

$  Modified:  Byoung  Yoon  Park  12-02-2011 

$  Experiment:  Hohler  &  Stilp  W  rod  vs  RHA 

$  Comments:  Validation  test  for  X-FEM  in  ALEGRA 

$===  VARIABLES==== 

$  {includeC'parameters.txt")} 

$  Time  Interval:  {T_INT  =  l.Oe-6}  s 

$  Plot  Interval:  {P_INT  =  T_INT}  s 

$  ===  execution  control  === 

$ - Initiation  and  Termination 

Title 

Hohler-Stilprha  {T_CAS},  v={V_ROD}  m/s,  Plate  Thick={T_PLT*m_cm}  cm,  CELL={S_CEL*m_mm}  mm 
units,  si 

READ  RESTART  DUMP  =  -1  $  Restart  at  the  latest  available  restart  dump 

termination  time,  {tt=160.0e-6}  $  Simulation  terminated  at  {tt*s_mus}  mus. 

TERMINATION  CPU  {tc=9.*h_s}  $  Run  will  be  terminated  within  CPU  time  {tc*s_h}  hr/64nodes(XFEM) 

Stermi nation  cycle,  1 

CHECK  SHUTDOWN  FILE  {csf=2}m  $  Check  shutdown  file  every  {csf}  minutes 

EMIT  RESTART,  TIME  INTERVAL  {T_INT}  $  Make  dmp  file  every  {T_lNT*s_mus}  mus. 

$ -  I/O  Control 

emit  plot,  time  i nterval={P_lNT}  $  s 

$emit  plot,  cycle  interval=10  from  time  7.e-5  to  15.e-5 

emit  hisplt,  cycle  interval=10 

$emit  hisplt,  time  i nterval={T_lNT}  $  s 

$emit  hisplt,  cycle  interval=l  from  time  7.e-5  to  15.e-5 

plot  variables 
no  underscores 
velocity 
density 
density, avg 
temperature 
temperature, avg 
pressure 
pressure , avg 
mat_mi n_coords 
mat_max_coo  rd  s 
end 

history  plot  variables 
no  underscores 
velocity 
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density 

temperature 

pressure 

end 


spy 

PlotTime(0.0,{P_lNT})  ; 

define  main() 

\{ 

pprintfC"  PLOT:  Cycle=%d,  Time=%e\n" .CYCLE, TIME) ; 

%  ListVariablesC)  ; 

Image ("Hohl ersti 1 pRha" .WHITE , BLACK) ; 

Label ("{T_CAS} ,  v={V_ROD}  m/s,  Plate  Thi ck={T_PLT*m_cm}  cm"); 

SMOOTH_SHADING  =  OFF; 

ColorMapRangeCO. ,20000. ,LIN_MAP) ; 

%  AutoScaleMapC'DENSiTY"  ,0.  ,20000.)  ; 

Drawcol orMap ("DENSITY" ,0.1, 0.3, 0.3, 0.8); 

Plot2D("DENSITY") ; 

Color(PURPLE) ; 

Draw2DMatContour ; 

Color(BLACK) ; 

DrawText(spri ntf ("t  =  %.2f  ~m~s",l.e6*TlME) ,0.4,0.05) ; 

Color(BLACK) ; 

Logoshadi ng(ON) ; 

LogoText(ON) ; 

LogoColors(BLACK, CYAN, BLACK) ; 

LogoPosi ti on (0 . 0 , 0 . 0 , 0 . 4 , 0 . 15) ; 

DrawSNL_Logo(l. 3) ; 

Endimage ; 

\} 

endspy 

$  ===  PHYSICS  INPUT=== 
solid  dynamics 
cyl  i  ndri  cal 

constant  volume  fraction  algorithm 
end 

xfem 

tangential  search  tolerance,  0.05 
maxenf  iterations,  100 
gap  rate  tolerance,  l.e-8 
mspair  2,1 
end 

domai n 

remap  method,  intersection 
pattern  reconstruction 
interface  order  201 
pattern  smoothing  =  off 
pattern  fragment  motion  =  off 
end 

time  step  scale,  {tss=.4}  $  {tss}  times  of  the  value  of  the  calculated  stable  time  step. 

$ - MESH - 

$BX  :{BX  =  1}  $  Number  of  Blocks  in  X 
$BY  :{BY  =  1}  $  Number  of  Blocks  in  Y 
$NX_TOT  :{NX_TOT  =  i nt ( (r_plt/s_CEL)/bx) } 

$NY_TOT  :{NY_TOT  =  i nt ( (L_MDL/S_CEL)/BY) } 

$NX  : {NX  =  NX_TOT/BX} 

$NY  : {NY  =  NY_TOT/BY} 

SNSLOTS  :{NSL0TS  =  BX*BY} 

$NUM_CELLS  :{NUM_CELLS  =  NX_TOT*NY_TOT} 

$R_DER  :{R_DER  =  4.*R_ROD}  $  Dense  element  radius  bound  (R-di r) 

$Z_DEU  :{Z_DEU  =  0.078}  $  Dense  element  upper  bound  (Z-di r)  (hole  depth  is  0.0741  at  3000  m/s) 

$Z_DEL  :{Z_DEL  =  0.000}  $  Dense  element  lower  bound  (Z-di r) 

$I_DER  :{I_DER  =  -100.}  $  Increase  Rate  (R-di r) ,  The  larger,  the  larger  increase  rate 

$I_DEU  :{I_DEU  =  -40.0}  $  Increase  Rate  (+Z-dir),  The  larger,  the  larger  increase  rate 

$I_DEL  :{I_DEL  =  70.0}  $  Increase  Rate  (-Z-dir),  The  larger,  the  smaller  increase  rate 

$F_DER  :{F_DER  =0.1}  $  Increase  Factor  (R-dir),  The  larger,  the  smaller  interval  at  end 

$F_DEU  :{F_DEU  =  0.9}  $  Increase  Factor  (+Z-di r) ,  The  larger,  the  smaller  interval  at  end 

$F_DEL  :{F_DEL  =  0.2}  $  Increase  Factor  (-Z-dir),  The  larger,  the  larger  interval  at  end 

$$$  T_PLT=160  mm,  l_DEU=-50,  S_CEL=0.15 

$$$  T_PLT=170  mm,  l_DEU=-45,  S_CEL=0.2 
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$$$  T_PLT=180  mm,  I_DEU=-40,  S_CEL=0.17 


mesh,  inline 
recti  1 i near 
bx  =  {BX} 
by  =  {BY} 
nx  =  {NX} 
ny  =  {NY} 

gmin  =  {C_MDL}  {F_MDL} 

gmax  =  {NX_TOT*BX*S_CEL}  {F_MDL+NY_TOT*BY*S_CEL} 
end 

set  assign 

sideset,ilo, 30 
si  deset, i hi ,10 
sideset, qlo,40 
sideset,]hi ,20 
end 

USER  DEFINED  ELEMENT  DENSITY,  X 
fi  el  d=l; 

if(coord  <  {R_DER})  {"{"}field  =  1;{"}"} 

ifCcoord  >=  {R_DER})  {"{''}field  =  1.  *exp({l_DER}*((coord)-{R_DER*F_DER}))  ;  {"}"} 

END 

USER  DEFINED  ELEMENT  DENSITY,  Y 
fi  el  d=l; 

if(coord  <  {Z_DEU})  {"{"}field  =  1;{"}"} 

ifCcoord  >=  {Z_DEU})  {"{"}field  =  1. *exp({l_DEU}*((coord)-{Z_DEU*F_DEU})) ; {"}"} 
ifCcoord  <=  {Z_DEL})  {"{''}field  =  l.*expC{l_DEL}*CCcoord)-{Z_DEU*F_DEL}))  ;{"}"} 
ifCcoord  <  {T_ROD})  {"{"}field  =  1;{"}"} 

END 

end 

$$ 

$  -  BOUNDARY  CONDITION  - 

no  displacement,  sideset  20, z  $  +Z 
no  displacement,  sideset  10, r  $  +R 
no  displacement,  sideset  30, r  $  -R  =  0 

block  1 

EULERIAN  MESH 
add  diatom  input 
end 

$  -  MATERIAL  INSERTION  - 

$ - 

$  material  1,  rod,  93%  WHA,  Tungsten 
$  material  2,  plate,  HzB  Steel 

$ - 

di atom 

translateCO. ,{H_ROD}) 
package  'ROD' 
material  1 
yvelocity  {V_ROD} 
insert  box 

pi  {C_MDL}  {H_ROD} 

p2  {R_ROD}  {T_ROD} 

endi 
endp 
endt 

translateCO. ,{F_PLT}) 
package  'PLATE' 
material  2 
insert  box 

pi  {C_MDL}  {F_PLT} 

p2  {R_PLT}  {B_PLT} 

endi 
endp 
endt 

enddiatom 

$  -  TRACERS  - 

tracer  points 

lagrangian  Tracer  1  x  {C_MDL}  y  {H_ROD}  $  Center  of  Rod  Head  --  Tracer  H 

lagrangian  tracer  2  x  {C_MDL}  y  {T_ROD}  $  Center  of  Rod  Tail 

lagrangian  tracer  3  x  {C_MDL}  y  {F_PLT}  $  Center  of  Plate  Front 

lagrangian  tracer  4  x  {R_PLT-S_CEL}  y  {F_PLT}  $  Edge  of  Plate  Front  --  Tracer  P 

lagrangian  tracer  5  x  {R_ROD}  y  {F_PLT}  $  Edge  of  Rod  on  Plate  Front 

end 
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plot,  exotracer 

file  =  'tracer.exo' 
end 

$  plot,  vtk  mesh 
$  file  =  'hohlerstilprha' 
$  end 

end 


$  ===  MATERIAL  === 
material  1  $  93%  WHA 

model  =  1 
end 

model  1  cth  elastic  plastic 
eos  model  =  100 

yield  model  =  10 

$  fracture  model  =  1000 
end 

model  10  johnson  cook  ep 

ajo  =  1.365e9  $  Hohler  &  Stilp 

bjo  =  0.1765e9 
cjo  =  0.016 
mjo  =  1.00 
njo  =  0.12 
tjo  =  3695. 
poisson  =  0.281 
end 

model  100,  keos  sesame 
feos  =  'sesame' 

neos  =  3550  $  TUNGSTEN 

sr  =  1.094  $  Scaled  for  17.760  g/cc 

clip  =  1.0 
end 


material  2  $  HzB  Steel 

model  =  2 
end 

model  2  cth  elastic  plastic 
eos  model  =  200 

yield  model  =  20 

end 

model  20  johnson  cook  ep 

ajo=0.810e9  $  Hohler  &  Stilp  HzB 

bjo=0.5095e9 
cj 0=0. 014 
mj 0=1. 030 
njo=0 . 260 
tj 0=1818 . 
poisson  =  0.299 
end 

model  200,  keos  sesame 
feos  =  'sesame' 
neos  =  2150  $  Iron 

sr  =  1.003  $  Scaled  for  7.85  g/cc 

clip  =  1.0 
end 

exi  t 
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Appendix  3:  Input  deck  for  5x12fixed 


[parameters.txt] 


Define  Parameters 

Test  Case: 

Cell  size: 

Model  Center  (X) : 

Radius  of  Plate  (X): 

Thickness  of  Plate  (Y) : 

Radius  of  Rod  (X): 

Lengh  of  Rod  (Y) : 

Velocity  of  Rod: 

Front  of  Plate  (Y) : 

Back  of  Plate  (Y)  : 

Head  of  Rod  (Y) : 

Tail  of  Rod  (Y) : 

Back  of  Model  (Y) : 

Front  of  1st  Graded  Area  (Y) : 
Radius  of  1st  Graded  Area  (X) : 


T_CAS  =  {T_CAS='5xl2fx0500'} 
S_CEL  =  {S_CEL=0.17*mm_m}  m 


C_MDL  =  {C_MDL= 
R_PLT  =  {R_PLT= 
T_PLT  =  {T_PLT= 
R_ROD  =  {R_ROD= 
L_ROD  =  {L_ROD= 
V_ROD  =  {V_ROD= 
F_PLT  =  {F_PLT= 
B_PLT  =  {B_PLT= 
H_ROD  =  {H_ROD= 
T_ROD  =  {T_ROD= 
B_MDL  =  {B_MDL= 
F_1GA  =  {F_1GA= 


0.00}  m 
50.0*mm_m}  m 
120.0*mm_m}  m 
2.5*mm_m}  m 
50.0*mm_m}  m 
500}  m/s 
0.00}  m 
F_PLT+T_PLT}  m 
F_PLT}  m 
H_ROD-L_ROD}  m 
R_PLT}  m 

78.*mm_m}  m  (hole  depth  is  0.0741  at  3000  m/s) 
4.*R_ROD}  m 


R_1GA  =  {R_lGA: 

S  Gap  between  Rod  tail  and  Fornt  of  Model  (Y) :  L_GAP  =  {L_GAP=  S_CEL}  m 
$  Front  of  Model  (Y) :  F_MDL  =  {F_MDL=  T_ROD-L_GAP}  m 

$  Length  of  Model  (Y) :  L_MDL  =  {L_MDL=  T_PLT+L_ROD+L_GAP}  m 

Note:  V_ROD  varied  from  500  m/s  to  3000  m/s  for  the  projectile  rod  velocities, 
[hohlerstilprha.inp] 

$  {i ncl ude("/home/bypark/common/uni ts . txt") } 

$  {_FORMAT="%.10g"} 


-  COMMENTS  -- 

Author: 

Modi fi ed : 
Experiment: 
Comments : 


Brian  Leavy 

Byoung  Yoon  Park  12-02-2011 
Hohler  &  Stilp  W  rod  vs  RHA 

Validation  test  for  X-FEM  in  ALEGRA 


$===  VARIABLES==== 

$  {includeC'parameters.txt")} 

$  Time  Interval:  {T_INT  =  l.Oe-6}  s 

$  Plot  Interval:  {P_INT  =  T_INT}  s 


$  ===  execution  control  === 

$ - Initiation  and  Termination 

Title 

Hohler-Stilprha  {T_CAS},  v={V_ROD}  m/s,  Plate  Thick={T_PLT*m_cm}  cm,  CELL={S_CEL*m_mm}  mm 
units,  si 

READ  RESTART  DUMP  =  -1  $  Restart  at  the  latest  available  restart  dump 

termination  time,  {tt=160.0e-6}  $  Simulation  terminated  at  {tt*s_mus}  mus. 

TERMINATION  CPU  {tc=5.*h_s}  $  Run  will  be  terminated  within  CPU  time  {tc*s_h}  hr/64nodes(XFEM) 

Stermi nation  cycle,  1 

CHECK  SHUTDOWN  FILE  {csf=2}m  $  Check  shutdown  file  every  {csf}  minutes 

EMIT  RESTART,  TIME  INTERVAL  {T_INT}  $  Make  dmp  file  every  {T_lNT*s_mus}  mus. 

$ -  I/O  Control 

emit  plot,  time  interval  ={P_INT}  $  s 

$emit  plot,  cycle  interval=10  from  time  7.e-5  to  15.e-5 

emit  hisplt,  cycle  interval=10 

$emit  hisplt,  time  i nterval={T_lNT}  $  s 

$emit  hisplt,  cycle  interval=l  from  time  7.e-5  to  15.e-5 

plot  variables 
no  underscores 
velocity 
density 
density, avg 
temperature 
temperature, avg 
pressure 
pressure , avg 
mat_mi n_coords 
mat_max_coo  rd  s 
end 


history  plot  variables 
no  underscores 
velocity 
density 
temperature 
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pressure 

end 


spy 

PlotTime(0.0,{P_lNT})  ; 

define  main() 

\{ 

pprintfC"  PLOT:  Cycle=%d,  Time=%e\n" .CYCLE, TIME) ; 

%  ListVariablesC)  ; 

Image ("Hohl ersti 1 pRha" .WHITE , BLACK) ; 

Label ("{T_CAS} ,  v={V_ROD}  m/s,  Plate  Thi ck={T_PLT*m_cm}  cm"); 

SMOOTH_SHADING  =  OFF; 

ColorMapRange(0. ,20000. ,LIN_MAP) ; 

%  AutoScaleMapC'DENSiTY"  ,0.  ,20000.)  ; 

Drawcol orMap ("DENSITY" ,0.1, 0.3, 0.3, 0.8); 

Plot2D("DENSITY") ; 

Color(PURPLE)  ; 

Draw2DMatContour ; 

Color(BLACK) ; 

DrawText(spri ntf ("t  =  %.2f  ~m~s" ,l.e6*TlME) ,0.4,0.05) ; 

Color(BLACK) ; 

Logoshadi ng(ON) ; 

LogoText(ON) ; 

LogoColors(BLACK, CYAN, BLACK) ; 

LogoPosi ti on (0 . 0 , 0 . 0 , 0 . 4 , 0 . 15) ; 

DrawSNL_Logo(l. 3) ; 

Endimage ; 

\} 

endspy 

$  ===  PHYSICS  INPUT=== 
solid  dynamics 
cyl  i  ndri  cal 

constant  volume  fraction  algorithm 
end 

xfem 

tangential  search  tolerance,  0.05 
maxenf  iterations,  100 
gap  rate  tolerance,  l.e-8 
mspair  2,1 
end 

domai n 

remap  method,  intersection 
pattern  reconstruction 
interface  order  201 
pattern  smoothing  =  off 
pattern  fragment  motion  =  off 
end 

time  step  scale,  {tss=.4}  $  {tss}  times  of  the  value  of  the  calculated  stable  time  step. 

$ - MESH - 

$BX  :{BX  =  1}  $  Number  of  Blocks  in  X 
$BY  :{BY  =  1}  $  Number  of  Blocks  in  Y 
$NX_TOT  :{NX_TOT  =  i nt ( (r_plt/s_CEL)/bx) } 

$NY_TOT  :{NY_TOT  =  i nt ( (L_MDL/S_CEL)/BY) } 

$NX  : {NX  =  NX_TOT/BX} 

$NY  : {NY  =  NY_TOT/BY} 

SNSLOTS  :{NSL0TS  =  BX*BY} 

$NUM_CELLS  :{NUM_CELLS  =  NX_TOT*NY_TOT} 

$R_DER  :{R_DER  =  4.*R_ROD}  $  Dense  element  radius  bound  (R-di r) 

$Z_DEU  :{Z_DEU  =  0.078}  $  Dense  element  upper  bound  (Z-di r)  (hole  depth  is  0.0741  at  3000  m/s) 

$Z_DEL  :{Z_DEL  =  0.000}  $  Dense  element  lower  bound  (Z-di r) 

$I_DER  :{I_DER  =  -100.}  $  Increase  Rate  (R-di r) ,  The  larger,  the  larger  increase  rate 

$I_DEU  :{I_DEU  =  -80.0}  $  Increase  Rate  (+Z-dir),  The  larger,  the  larger  increase  rate 

$I_DEL  :{I_DEL  =  70.0}  $  Increase  Rate  (-Z-dir),  The  larger,  the  smaller  increase  rate 

$F_DER  :{F_DER  =0.1}  $  Increase  Factor  (R-di r) ,  The  larger,  the  smaller  interval  at  end 

$F_DEU  :{F_DEU  =0.7}  $  Increase  Factor  (+Z-di r) ,  The  larger,  the  smaller  interval  at  end 

$F_DEL  :{F_DEL  =  0.3}  $  Increase  Factor  (-Z-dir),  The  larger,  the  larger  interval  at  end 

$$$  T_PLT=160  mm,  R_PLT=  40,  l_DEU=-50,  S_CEL=0.15 
$$$  T_PLT=170  mm,  R_PLT=  40,  I_DEU=-45,  S_CEL=0.2 
$$$  T_PLT=180  mm,  R_PLT=  40,  I_DEU=-40,  S_CEL=0.17 
$$$  T_PLT=120  mm,  R_PLT=  50,  I_DEU=-80,  S_CEL=0.17 
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mesh,  inline 
recti  1 i near 
bx  =  {EX} 
by  =  {BY} 
nx  =  {NX} 
ny  =  {NY} 

gmin  =  {C_MDL}  {F_MDL} 

gmax  =  {NX_TOT*BX*S_CEL}  {F_MDL+NY_TOT*BY*S_CEL} 
end 

set  assign 

sideset.ilo, 30 
si  deset, i hi ,10 
sideset, jlo,40 
sideset, jhi ,20 
end 

USER  DEFINED  ELEMENT  DENSITY,  X 
fi  el  d=l; 

if(coord  <  {R_DER})  {"{"}field  =  1;{"}"} 

ifCcoord  >=  {R_DER})  {"{"{field  =  1. *exp({l_DER}*C(coord)-{R_DER*F_DER})) ; {"}"} 

END 

USER  DEFINED  ELEMENT  DENSITY,  Y 
fi  el  d=l; 

if(coord  <  {Z_DEU})  {"{"{field  =  1;{"}"} 

ifCcoord  >=  {Z_DEU})  {"{"{field  =  l.*exp({l_DEU}*((coord)-{Z_DEU*F_DEU})) ;{"}"} 
ifCcoord  <=  {Z_DEL})  {"{"{field  =  1. *expC{l_DEL}*CCcoord)-{Z_DEU*F_DEL})) ; {"}"} 

ifCcoord  <  {t_rod}3  {"{"{field  =  1;{"{"} 

END 

end 

$$ 

$  -  BOUNDARY  CONDITION  - 

no  displacement,  sideset  20, z  $  +Z 
no  displacement,  sideset  10, r  $  +R 
no  displacement,  sideset  30, r  $  -R  =  0 

block  1 

EULERIAN  MESH 
add  diatom  input 
end 

$  -  MATERIAL  INSERTION  - 

$ - 

$  material  1,  rod,  93%  WHA,  Tungsten 
$  material  2,  plate,  HzB  Steel 

$ - 

di atom 

translateCO. ,{H_ROD{) 
package  'ROD' 
material  1 
yvelocity  {V_ROD{ 
insert  box 

pi  {C_MDL{  {H_ROD{ 

p2  {R_ROD}  {T_ROD} 

endi 
endp 
endt 

translateCO. ,{F_PLT{) 
package  'PLATE' 
material  2 
insert  box 

pi  {C_MDL{  {F_PLT{ 

p2  {R_PLT}  {B_PLT} 

endi 
endp 
endt 

enddiatom 

$  -  TRACERS  - 

tracer  points 

lagrangian  Tracer  1  x  {C_MDL{  y  {H_ROD{  $  Center  of  Rod  Head  --  Tracer  H 

lagrangian  tracer  2  x  {C_MDL}  y  {T_ROD}  $  Center  of  Rod  Tail 

lagrangian  tracer  3  x  {C_MDL}  y  {F_PLT}  $  Center  of  Plate  Front 

lagrangian  tracer  4  x  {R_PLT-S_CEL{  y  {F_PLT}  $  Edge  of  Plate  Front  --  Tracer  P 

lagrangian  tracer  5  x  {R_ROD{  y  {F_PLT}  $  Edge  of  Rod  on  Plate  Front 

end 

plot,  exotracer 
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file  =  'tracer.exo' 
end 

$  plot,  vtk  mesh 
$  file  =  'hohlerstilprha' 
$  end 

end 


$  ===  MATERIAL  === 
material  1  $  93%  WHA 

model  =  1 
end 

model  1  cth  elastic  plastic 
eos  model  =  100 

yield  model  =  10 

$  fracture  model  =  1000 
end 

model  10  johnson  cook  ep 

ajo  =  1.365e9  $  Hohler  &  Stilp 

bjo  =  0.1765e9 
cjo  =  0.016 
mjo  =  1.00 
njo  =  0.12 
tjo  =  3695. 
poisson  =  0.281 
end 

model  100,  keos  sesame 
feos  =  'sesame' 

neos  =  3550  $  TUNGSTEN 

sr  =  1.094  $  Scaled  for  17.760  g/cc 

clip  =  1.0 
end 


material  2  $  HzB  Steel 

model  =  2 
end 

model  2  cth  elastic  plastic 
eos  model  =  200 

yield  model  =  20 

end 

model  20  johnson  cook  ep 

ajo=0.810e9  $  Hohler  &  Stilp  HzB 

bjo=0.5095e9 
cj 0=0. 014 
mj 0=1. 030 
njo=0 . 260 
tj 0=1818 . 
poisson  =  0.299 
end 

model  200,  keos  sesame 
feos  =  'sesame' 
neos  =  2150  $  Iron 

sr  =  1.003  $  Scaled  for  7.85  g/cc 

clip  =  1.0 
end 

exi  t 
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